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1  -  INTRODUCTION 


The  code  is  named  CANTATA  -  CANard/TAil  Transonic  Aerodynamics.  It 
operates  on  the  CRAY  X-MP  computer  at  NASA  Ames  Research  Center.  Part  of  the 
computer  time  for  the  code  development  and  evaluation  was  provided  through  the 
Applied  Computational  Fluids  Branch  at  Ames.  The  code  is  stored  on  the  Ames 
CRAY  in  UPDATE  format. 

The  CANTATA  code  was  developed  from  the  PANDORA  wing-body-canard  code 
(Ref  1).  The  PANDORA  code  was  developed  from  an  earlier  version  (1978)  of  the 
Boppe  wing-body  code  (Ref  2).  Some  parts  of  the  three  codes  are  essentially 
identical.  Unnecessary  details  of  the  code  which  are  the  same  as  described  in 
Ref  1,  2  and  4  are  not  repeated  here.  Those  references  should  be  consulted  if 
the  particular  code  details  are  of  interest. 

The  major  improvements  incorporated  in  the  CANTATA  code  are  arbitrary 
canard  (or  tail)  placement,  an  Alternating  Direction  Implicit  (ADI)  solution 
and  a  "floating”  wake  analysis  capability.  In  addition,  the  CANTATA  code 
takes  advantage  of  the  vector  processing  capability  of  the  CRAY  computer. 


1  Aidala,  P. ,  "Numerical  Aircraft  Design  Using  3-D  Transonic  Analysis  with 
Optimization,"  AFWAL-TR-8 1-3091 ,  August  1981. 

2  Boppe,  C.  W.,  "Transonic  Flowfield  Analysis  for  Wing-Fuselage  Combina¬ 
tions,"  NASA  Contractor  Report  3243,  May  1980. 


2  -  WING-BODY-CANARD  TRANSONIC  ANALYSIS 


The  basic  computational  method  employed  in  the  transonic  analysis  is  that 
of  Boppe  (Ref  2).  The  significant  new  developments  are  described  in  detail 
here.  The  common  parts  (flow  equation,  embedded  grid  interfacing,  body  model¬ 
ing,  viscous  effects)  are  described  briefly.  Additional  details  can  be  found 
in  Ref  1  and  2. 

2.1  FLOW  EQUATION 

The  flow  equation  used  in  the  analysis  is  an  "extended"  small-disturbance 
equation: 

(i-MZ  -  (v+i  )m24>x  -  m2*x>*xx  '  2h2<py 

+  (i-(y-Dm2<|>x)<|>yy  +  <|>2Z  -  0 

The  additional  terms  have  been  added  to  better  capture  swept  shock  waves  and 
more  accurately  determine  the  critical  velocity.  Empirical  modifications  and 
similarity  variables  are  not  employed.  Pressure  coefficients  on  wing  surfaces 
are  computed  using  the  complete  isentropic  formula.  To  simplify  velocity  com¬ 
putations  on  the  non-planar  body  surface,  a  simplified  equation  is  used: 

C  -  -  (2<J)  +  (i-m2)4>2) 

p  X  X  (2) 

The  computational  space  is  filled  with  a  relatively  crude  Cartesian  mesh. 
Instead  of  adopting  a  formal  far-field  solution  for  the  grid  outer  boundaries, 
the  original  X,  Y,  Z  region  is  stretched  to  i,  r|»  C>  a  region  in  which  the 
boundaries  correspond  to  infinity.  The  flow  field  potential  is  set  to  zero  on 
all  bounding  planes  except  the  downstream  plane,  at  which  all  X-derivatives 
are  set  to  zero, 
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The  following  conditions  are  enforced  at  the  symmetry  plane. 

4>y  =  o 

*XY  -  ° 

2.2  GLOBAI  CRUDE  GRID 

The  key  item  to  enable  wing-canard  analysis  without  undue  computer  stor¬ 
age  requirements  was  the  development  of  suitable  grid  point  distributions. 

The  global  crude  mesh  uses  65,  26  and  31  planes  in  the  streamwise,  spanwise 
and  vertical  directions,  respectively.  This  results  in  52,390  mesh  points. 

In  order  to  distribute  the  crude  grid  streamwise  mesh  planes  as  effec¬ 
tively  as  possible,  separate  grids  are  used  for  wing  and  wing-canard  analysis. 
Without  a  canard,  a  symmetric  streamwise  grid  is  generated  with  the  function 

X  -  Bq£  (4a) 

X  -  TAN  (3  (£  -  £j)  )  +  B2TAN  (|  (|  -  (4b) 

where  Bq  is  chosen  to  place  80%  of  the  grid  points  between  ±X^.  The  value  of 
Xj  is  chosen  to  center  the  grid  between  the  most  forward  and  aft  points  of  the 
exposed  wing.  The  constant  B1  is  used  to  match  the  grid  spacing  (dx/d£)  at 
Xj ,  and  B2  is  used  to  provide  additional  stretching  to  the  far  field. 

The  length  scale  for  the  streamwise  grid  is  the  chord  at  midspan  of  the 
wing.  The  transformation  places  the  first  and  last  points  at  infinity.  The 
second  and  f irst-to-last  points  are  the  effective  numerical  boundaries,  and 
are  at  8.6  mid-chords  from  the  mid-point  of  the  chord  at  mid-span.  The 
far-field  boundary  is  not  swept  or  tapered. 

When  a  canard  is  present,  the  grid  of  Equation  (4  a,b)  usually  results  in 
too  few  crude  grid  planes  intercepting  the  canard  planform.  The  streamwise 
grid  transformation  used  for  wing-canard  combinations  is: 

X  -  AQ  +  AjC  +  A212  +  A3|3  +  aJ?  +  A5|5  +  AsC/(l-f) 

l  -  C/4* 


(3a) 

(3b) 


(5a) 

(5b) 


This  transforms  the  finite  domain  ~£  <_ £  £  4^  to  the  infinite  region 

-oo  <  x  <  “.  Constant  is  empirical,  controlling  the  rate  of  stretching  near 
infinity.  The  points  at  £,  ~  ±  1  are  made  the  points  of  maximum  density  of 
mesh  planes  by  specifying  the  second  derivative  of  X  to  be  zero.  The  value  of 
X  and  the  first  derivative  of  X  at  |  =  ±  1  are  also  specified  to  determine  the 
coefficients  Aq  to  A^.  The  value  of  4^  is  adjusted  iteratively  to  place  ap¬ 
proximately  607.  of  the  total  mesh  planes  between  the  most  forward  and  the  most 
aft  points  on  the  wing-canard  combination. 

For  an  aft  swept  wing-canard  combination,  the  wing  tip  and  canard  tip 
determine  the  transformation.  The  mid-points  of  the  two  tip  chords  are  the 
values  of  X  for  |  =  i  1.  The  first  derivative  of  X  is  set  to  result  in  a 
nominal  six  mesh  planes  intercepting  each  of  the  tip  chords.  An  example  of 
the  resulting  wing-canard  grid  transformation  and  the  corresponding  physical 
mesh  plane  distribution  is  shown  in  Figure  1.  The  mesh  generation  has  been 
applied  to  several  aft-swept  configurations,  several  forward-swept  configura¬ 
tions  and  to  several  arbitrary  wing-canard  parametric  variations.  Good  re¬ 
sults  were  observed  in  all  cases.  For  forward  swept  wings,  the  canard  tip  and 
wing  tip  may  be  at  the  same  streamwise  location.  In  this  case,  the  two 


Downstream  infinity  Physical  space  grid 

Figure  1.  Wing-canard  crude  mesh  transformation. 


streamwise  points  used  for  ^  =  t  1  in  the  transformation  are  the  wing/canard 
tip  location  and  near  the  trailing  edge  of  the  centerline  of  the  wing. 

The  two  arbitrary  factors,  Ag  and  4^  do  not  allow  an  analytic  expression 
for  the  distance  to  the  far  field.  In  addition,  the  local  spacing  is  adjusted 
for  the  canard  or  wing  chord,  so  that  the  grid  is  not  symmetrical.  As  ex¬ 
plained  in  the  above  discussion,  the  grid  scales  with  the  two  tip  chords  and 
their  streamwise  separation.  This  is  independent  of  the  spanwise  location  of 
the  tips. 

The  spanwise  grid  transformation  is  much  simpler: 

Y  =  (l/A2)TANH_1(n)  +  +  C3n3  (6) 

This  transforms  the  domain  0<ri<lto0<Y<®.  Constant  A.^  is  determined 
to  place  18  spanwise  mesh  planes  between  the  centerline  and  wing  tip,  with  the 
wing  tip  falling  midway  between  computational  planes.  Constants  Cj  and  C^  are 
used  to  perturb  the  grid  so  that  the  canard  tip  falls  midway  between  two  mesh 
planes  and  the  wing  tip  location  is  not  changed.  Constants  C^  and  C^  are  set 
to  zero  when  a  canard  is  not  present.  For  a  wing-only  grid,  the  spanwise 
boundary  (next~to~last  point)  is  at  2.2  semispans.  When  a  canard  is  present, 
the  change  in  the  spanwise  boundary  is  usually  insignificant. 

The  vertical  grid  transformation  combines  a  constant  spacing  interval 
near  the  lifting  surfaces  with  a  tangent  function  stretching  to  the  far  field. 
The  computational  domain  -1  <  £  £  1  is  transformed  to  the  physical  domain 
through  the  relation 


MID  +  B0^ 

^LO  —  ^  —  ^up 

(7a) 

MID  +  BO^UP  +  B1^"^UP^  +  B3TAN^B4^_^up^ 

t  _  £UP 

(7b) 

MID  +  B0^L0  +  B1^-^L0^  +  b3tan(B4(£~£lo^ 

<T* 

1  A 

(Tt 

o 

(7c) 

B0  -  B3B4 

(7d) 

bo^up  +  Azbndry  "  ^1-^”^UP^B0 

TAN(B4(1-A£-Cup)  -  B4(1-AHUP) 

tit. 

(7e) 

BO^LO  +  AZBNDRY  ~  0_AC+Clo^bo 

tan(b4u-aC+Cl0>  -  BA(1-ACHL0) 

o 

*-/» 

V| 

(70 

(7f) 


The  coefficient  Bq  is  initially  chosen  for  a  grid  spacing  of  0.088  mid-chords. 

It  is  then  adjusted  to  align  grid  planes  with  the  two  lifting  surfaces,  with 

the  center  of  the  grid  midway  between  the  surfaces  (at  Z  ).  The  constant 

MID 

spacing  of  Equation  (7a)  is  extend  two  grid  planes  above  and  below  the  sur¬ 
faces  to  and  respectively.  The  stretchings  of  Equations  (7b)  and 

(7c)  place  the  next  to  last  points  a  distance  AZ  ,  which  is  eight  mean 

dNDKY 

aerodynamic  chords,  from  Z  . 

MID 

An  input  parameter  (ITWAKE)  is  used  to  turn  on  the  free  wake  calculation. 
Details  of  the  calculation  are  given  in  Appendix  A.  A  wake  shape  calculation 
is  performed  every  ITWAKE  iterations  of  the  potential  solution.  After  each 
wake  shape  calculation,  the  Z  grid  is  sheared  to  follow  the  vortex  sheet  that 
leaves  the  trailing  edges  of  the  lifting  surfaces.  The  wake  shape  of  the 
shorter  semispan  is  extended  spanwise  parallel  to  the  larger  semispan.  The 
larger  semi  span  is  extended  spanwise  at  a  constant  value  of  Z  equal  to  the 
most  outboard  value  on  the  wake.  The  Z  grid  is  then  reconstructed,  using  the 
same  stretching  formula  described  previously.  The  resulting  Z  grid  is  sheared 
in  X  and  Y,  and  does  not  permit  two  wakes  to  merge. 

2.3  EMBEDDED  FINE  GRID  SYSTEM 

Individual  fine  grid  arrays  are  constructed  for  the  wing  and  canard.  The 
code  always  designates  the  forward  lifting  surface  as  the  "canard"  and  the  aft 
surface  as  the  "wing."  Thus,  the  wing  of  a  wing-tail  configuration  is  con¬ 
sidered  the  "canard"  and  the  tail  is  considered  the  "wing."  These  secondary 
mesh  systems  serve  two  purposes.  First,  detailed  computations  are  performed 
only  in  a  region  very  close  to  the  surface  where  gradients  are  large  and  de¬ 
tails  are  important.  The  resulting  numerical  efficiency  permits  a  very  dense 
computational  mesh,  a  benefit  in  both  the  resolution  of  shock  waves  and  the 
calculation  of  configuration  forces  and  moments.  Second,  the  embedded  mesh 
systems  are  independent  and  optimized  for  a  particular  geometric  component. 

The  system  is  not  constrained  by  a  single  geometry-fitting  transformation. 

Fine  grid  arrays  are  set  up  at  each  position  where  a  crude  spanwise  mesh 
plane  cuts  the  wing  and  canard  surface.  With  18  spanwise  crude  planes  between 
the  centerline  and  wing  tip,  an  isolated  wing  analysis  would  have  18  fine  wing 
grids.  (A  fuselage  would  reduce  the  number  of  fine  wing  grids  by  the  number 


of  mesh  planes  within  the  computational  wing  root  junction.)  The  number  of 
canard  fine  grids  is  proportional  to  the  extent  of  the  canard  semispan.  If 
the  canard  semispan  were  half  that  of  the  wing,  then  the  canard  would  have 
approximately  half  the  number  of  wing  fine  grids.  This  would  increase  the 
computation  time  for  the  fine  grid  solution  by  50%. 

The  local  section  leading  edge  is  placed  midway  between  mesh  points  and 
the  trailing  edge  is  placed  at  a  mesh  point.  The  streamwise  mesh  spacing  is 
scaled  by  the  local  chord.  Thus,  each  section  has  the  same  number  of  grid 
points  along  the  chord  (72).  The  nominal  boundaries  of  the  fine  grids  are  at 
20%  of  the  local  chord  in  front  of  each  leading  edge  and  40%  behind  each 
trailing  edge.  The  total  number  of  fine  grid  points  will  vary  for  different 
canard/wing  semispan  ratios.  The  dimensions  of  the  code  will  allow  101,250 
fine  grid  points.  This  corresponds  to  30  fine  grids  on  the  wing  and  canard, 
with  each  having  135  streamwise  points  and  25  vertical  points. 

The  starting  fine  grid  system  is  illustrated  in  Figure  2.  As  described 
above,  the  fine  grids  are  sheared  and  tapered  to  conform  to  the  wing  and  canard 
planforms.  This  represents  a  transformation  function  between  the  fine  grid 
computational  domain  £  and  the  physical  domain  X: 

X  '  5  <XGTE  *  XGLE>  *  XGLE  <«> 
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Figure  2.  Embedded  fine  grid  system. 
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The  spanwise  (Y)  transformation  is  the  same  as  that  for  the  crude  grid.  The 
fine  Z  grid  is  aligned  with  the  crude  Z  grid  along  any  surfaces  and  wakes.  In 
addition,  the  top  and  bottom  of  the  fine  grid  is  aligned  with  the  appropriate 
crude  grid  plane.  Linear  interpolation  is  used  in  the  streamwise  direction. 
The  fine  grid  spacing  is  constant  between  surfaces.  The  constant  spacing  is 
continued  two  points  above  and  below  the  surfaces.  Then  the  grid  is  stretched 
with  a  parabolic  function  to  the  top  and  bottom  boundaries.  The  nominal  loca¬ 
tions  of  the  boundaries  are  five  crude  grid  planes  above  the  surface,  and 
three  crude  grid  planes  below  the  surface. 


2.4  SOLUTION  TECHNIQUE 

The  crude  grid  solution  uses  Successive  Line  Over-Relaxation  (SLOR).  As 
the  vertical  line  relaxation  marches  downstream,  each  streamwise  location 
(streamwise  index)  is  tested  to  see  if  a  surface  or  wake  will  be  intercepted. 
The  test  is  done  for  both  a  "high"  and  "low"  surface  or  wake.  The  vertical 
line  segments  are  subdivided  by  each  surface  or  wake.  Thus,  when  two  surfaces 
(or  wakes)  are  present  at  a  streamwise  index,  three  separate  tridiagonal 
inversions  occur.  When  the  two  surfaces  are  in-plane,  the  "low"  surface  is 
the  forward  one,  and  the  "high"  surface  is  the  aft  one,  with  the  appropriate 
logic  testing  for  the  in-plane  arrangement.  The  crude  grid  solution  starts 
at  the  spanwise  boundary  and  then  marches  toward  the  centerline.  The  fine 
grid  solution  update  is  developed  separately  on  the  wing  and  canard.  The 
solution  is  first  done  on  the  wing  and  then  on  the  canard.  When  the  fine  grid 
solution  on  the  canard  has  been  updated,  one  iteration  of  the  fine  solution  is 
completed . 


The  fine  grid  solution  of  Equation  (1)  employs  an  approximate  factoriza¬ 
tion  algorithm  that  was  suggested  by  Professor  Antony  Jameson  of  Princeton 
University  as  a  consultant  to  this  study.  The  factorization  can  be  written 
as : 


(a  ♦  BjTfi*)  [(a  +  B26y)  (a  +  B^)  +  E)(1]60n  -  wa2L<t>r 


4>n*1  =  0n  +  6<J>n 
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The  operator  L  produces  the  residual  of  Equation  (1).  The  coefficient  T  is 
the  nonlinear  coefficient  of  in  Equation  (1).  The  6*  operator  is  used  at 
elliptic  (subsonic)  points  in  the  flowfield,  and  the  6X  operator  is  used  at 
hyperbolic  (supersonic)  points. 

The  solution  of  this  factorization  is  developed  in  three  steps: 


Step 

1. 

(a  + 

BlT6x)6<l>l  “  m2l^n 

(10a) 

Step 

2. 

(o  + 

B2<5y)<5<J>2  “  -  Exl64>n 

(10b) 
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3. 

(a  + 

B 36|)6<Dn  -  fidij 

(10c) 

Each  step  involves  only  tridiagonal  inversions,  and  can  be  vectorized  by  in¬ 
verting  several  matrices  in  parallel.  Within  the  separate  wing  or  canard  fine 
grid  system,  Step  1  is  solved  for  all  points  in  the  flowfield.  Then  Steps  2 
and  3  are  marched  downstream. 

The  embedded  fine  grid  interface  with  the  crude  grid  is  accomplished  by 
alternately  updating  the  crude  grid  and  fine  grid  solutions.  Potential  values 
at  interior  points  (i.e.,  on  a  lifting  surface)  of  the  crude  grid  are  fixed  by 
interpolating  the  most  recent  fine  grid  solution.  The  potential  values  at  the 
perimeter  of  the  fine  grids  are  fixed  by  interpolating  the  most  recent  crude 
grid  solution.  In  order  to  speed  up  the  overall  solution  convergence,  the 
fine  grid  solution  is  not  calculated  until  the  crude  grid  solution  has 

established  the  "coarse"  characteristics  of  the  flow  (approximately  100  itera¬ 
tions).  Then  the  crude/fine  interaction  is  begun  until  both  grids  are  sat¬ 
isfactorily  converged  (typically  an  additional  105  cycles).  Three-point 
LaGrange  interpolation  of  the  crude  grid  potential  values  is  used  to  initial¬ 
ize  and  update  the  fine  grid. 

One  cycle  of  the  crude/fine  grid  interaction  consists  of  two  steps  (see 
Figure  3): 
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Figure  3.  Fine/crude  grid  interface  (from  Ref  2). 


The  embedded  wing  grid  is  swept  holding  fine  grid  perimeter 
points  fixed  as  an  outer  boundary.  Conventional  Neumann  bound¬ 
ary  conditions  (4>n)  are  imposed  at  fine  grid  section  boundary 
points  forming  an  inner  boundary. 

The  crude  grid  section  boundary  points  are  computed  using  the 
potentials  at  the  fine  grid  section  boundary  points  (linear  in¬ 
terpolation).  These  crude  potentials  ( <t» )  are  held  fixed  for  the 
global  crude  grid  sweep  forming  an  array  of  Dirichlet  inner 
boundary  conditions.  Infinity  boundary  conditions  at  the  limits 
of  the  crude  computational  space  form  the  outer  boundary.  At 
the  end  of  the  crude  grid  sweep,  crude  grid  potentials  are  used 
to  update  the  fine  grid  perimeter  points. 


One  solution  iteration  of  only  the  crude  grid  takes  0.35  sec  of  CRAY  X-MP 
CPU  time  (6.7x10  ^  CPU  sec  per  point).  One  crude/fine  cycle  takes  1.4  sec  CPU 
time  for  23  wing  and  canard  fine  grids.  Each  fine  grid  has  135  streamwise 
points  and  25  vertical  points.  Thus,  one  crude/fine  solution  cycle  solves 
130,015  points  total  (1.9x10  ^  CPU  sec  per  point).  For  different  wing-canard 
combinations ,  the  crude/fine  CPU  time  is  essentially  proportional  to  the  num¬ 
ber  of  wing  and  canard  fine  grid  arrays. 

The  development  of  the  solution  is  monitored  by  calculating  the  maximum 
update  to  the  flow  field  potential  (CCMAX  in  the  printout)  and  its  position  in 
the  three-dimensional  flow  field.  The  magnitude  of  the  potential  update  is 
inversely  proportional  to  the  acceleration  parameter  a  in  Equation  (10).  The 
parameter  is  cycled  through  eight  values  during  the  solution.  In  addition,  an 
overall  scaling  of  the  a  parameters  occurs  during  the  solution  based  on  the 
change  in  the  maximum  residual.  The  initial  values  of  the  parameter  are  set 
to  relatively  large  values  (i.e.,  slower  convergence)  in  order  to  provide  more 
reliable  results  for  difficult  cases.  The  automatic  scaling  will  later  de¬ 
crease,  or  increase,  the  a  parameters  as  appropriate.  However,  the  scaling 
occurs  only  after  each  set  of  eight  values,  which  corresponds  to  eight  solu¬ 
tion  iterations.  This  results  in  a  relatively  slow  decrease  of  the  or  param¬ 
eters  when  it  may  occur.  Consequently,  "easy"  solutions  are  penalized  with  a 
slower  rate  of  convergence  than  would  be  optimal. 

The  solution  is  essentially  converged  when  the  circulation  (CIR  in  the 
printout)  shows  no  significant  change  in  about  24  solution  cycles.  It  is  pos¬ 
sible  that  a  non-decreasing  CCMAX  value  is  due  to  "noise"  in  the  crude/fine 
interface  that  does  not  decay.  The  position  of  CCMAX  is  useful  in  pinpointing 
the  problem  area  when  divergence  occurs.  Problem  areas  may  develop  when  the 
flow  conditions  are  extreme  or  when  the  geometric  representation  is  in  error. 
The  solution  will  terminate  itself  if  divergence  is  detected. 

Three  other  parameters  that  indicate  the  solution  convergence  appear  in 
the  code  output  as  CCAV,  RSD  and  RSDAV.  The  CCAV  parameter  is  the  average  of 
the  absolute  value  of  the  flow  potential  corrections  at  every  point.  It  is  a 
more  reliable  indication  of  the  solution  convergence  than  CCMAX,  and  is  used 
to  satisfy  the  convergence  test  input  option.  The  RSD  parameter  is  the  larg¬ 
est  residual  in  the  solution.  The  location  of  the  RSD  value  is  also  printed. 


The  RSDAV  parameter  is  the  average  of  the  absolute  value  of  the  residual  at 
every  point.  The  vertical  index  for  the  location  of  CCMAX  and  RSD  are  nega¬ 
tive  when  the  point  is  located  at  the  lower  side  of  the  wing  or  canard  plane. 
The  jump  in  potential  at  these  locations  is  handled  by  an  extra  array  of  po¬ 
tential  values  that  have  their  own  index  system.  As  discussed  for  CCMAX,  the 
values  of  CCAV,  RSD  and  RSDAV  will  vary  due  to  changes  in  the  acceleration 
parameter  a. 

No  precise  non-dimensional  parameter  is  available  to  compare  the  solution 
convergence  for  different  grids.  The  CCMAX  and  CCAV  values  are  scaled  by  the 
wing  average  cord.  This  results  in  the  same  values  for  CCMAX  and  CCAV  if  the 
scale  of  the  input  geometry  is  changed  (e.g.,  full  scale  or  model  scale). 
Changes  to  the  planform  (e.g.,  aspect  ratio  or  taper)  would  result  in  CCMAX 
and  CCAV  values  that  are  not  directly  comparable.  As  mentioned  above,  the 
reduction  in  CCMAX  and  CCAV  relative  to  the  starting  value  should  be  used  to 
evaluate  the  solution  convergence. 

2.5  TRANSFORMED  QUATION  AND  FINITE  DIFFERENCE  APPROXIMATIONS 

Due  to  the  grid  stretching  and  shearing,  the  derivatives  in  terms  of  the 
physical  domain  coorcinaL.es  (X,Y,Z)  are  calculated  by  the  following  formulas 
used  with  the  compuc  tional  domain  variables  (t.fj.C): 


<t>x  <t>4ex  +  04dx 

(11a) 

*Y  “  <^Y  +  V*  + 

(lib) 

*Z  =  Z 

(11c) 

^xx  '  ^xx  +  ^xx  +  ^4|^x  +  2(|)|C^x^x  +  ^^x 
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^XY  =  \  ^XY  +  ^XY  +  ^X  ^Y  +  ^X  nY 

(Ilf) 
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In  the  crude  grid  the  X  grid  is  independent  of  Y,  so  that  the  £y. 

4xy.  and 

terms  are  not  needed.  They  are  not  coded  for  the  crude  grid.  The  £y, 

L  ,  C  ,  and  CVv  terms  are  zero  when  there  is  no  wake  deflection  (the  Z  grid 
AA  ai  i  y 


is  then  independent  of  X  and  Y).  The  shearing  terms  for  the  X  and  Z  grid  are 
calculated  with  second  order  differences  and  implicit  differentiation 
formulas . 


The  finite  difference  approximations  to  the  first  derivatives  in  the 
above  expressions  are  always  second-order,  e.g.: 


<t>t  =  (<!>k+1  "  <t>k_1)/(2AO  (12) 

When  the  flow  is  subsonic  all  of  the  second  derivatives  are  second  order, 
e.g.  : 

♦«k  -  <Vi  -  2\  *  Vi>"*2  <13> 

The  only  change  when  the  flow  is  supersonic  is  "upwind  differencing"  of  the 
streainwise  second  derivative: 

"  ( 2<t> j  ”  <1^  ~  2,l)j_1  +  <t>I_2)/AC2  O*) 

where  the  +  superscript  indicates  new  values  in  the  crude  grid  SLOR  solution. 
The  upwind  differencing  is  used  only  for  the  <|>||  contribution  to  (j)^.  The 
term  for  4>yy  is  always  central  differenced. 

The  crude  grid  overrelaxation  of  the  solution  is  incorporated  at  subsonic 
points  as: 

hi  "  <Vi  -  s  -  2<4>  *i  *  <>5> 

”1 

with  iu  usually  1.6.  The  spanwise  differences  make  use  of  the  new  potential 
values  at  J  +  1: 


*4n  ^I-l.J+l  ^I+l.J-1  +  ^1+1, J+l)/(4A^An) 

$nn  “  (4>I, J-1  '  *i,j  "  0I,J  +  ♦z,j+i)/An2 


(16) 

(17) 


The  vertical  line  relaxation  implicitly  uses  all  new  potential  values  for 

Boundary  conditions  are  imposed  by  setting  the  value  of  a  potential  or 
its  first  derivative  at  a  field  point  which  represents  the  conf iguration  sur¬ 
face.  The  wing  and  canard  are  approximated  by  planar  surfaces.  As  discussed 
previously,  the  wakes  can  be  either  fixed  plar ar  surfaces  or  "floating"  sur¬ 
faces  updated  periodically  during  the  solution.  The  body  or  fuselage  is  rep- 


resented  by  a  fixed  cross-sectional  surface  extending  from  upstream  to  down¬ 
stream  infinity.  Corrections  are  applied  to  the  body  boundary  conditions  for 
the  simulation  of  finite  length  bodies.  Modeling  is  sufficiently  flexible  to 
permit  the  treatment  of  wings  at  varying  height  relative  to  the  body  (high- 
low  in  i  d  wing). 

Special  difference  approximations  are  required  at  boundary  points.  The 
lifting  surfaces  and  wakes  are  represented  numerically  by  a  grid  surface  of 
double  valued  potentials.  For  a  wing  surface  defined  by 

Z  =  F ( X , Y )  (18) 
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the  wing  flow  tangency  condition  is  approximated  by 

0Z  (X , Y, 0)  =  Fx  ”  a  +  6x  (19 

.)V 

where  the  slope  of  the  boundary  layer  displacement  thickness,  o^,  is  added 
only  for  the  i nv i sc id/vi scous  interaction  mode  of  operation. 

The  wing  upper  and  lower  surface  boundary  conditions  enter  the  solution 
formulation  by  way  of  the  <h  term  in  Equation  (1).  At  a  surface,  this  term 


can  be  written 
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where  the  I  and  J  subscripts  have  been  dropped  for  convenience.  By  incorpo¬ 
rating  the  following  relation 
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the  wing  boundary  condition  on  the  upper  surface  becomes 
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Similarly,  the  wing  lower  surface  boundary  condition  becomes 


At  the  end  of  each  sweep  of  the  flow  field,  the  Kutta  condition  is  en¬ 
forced  by  calculating  the  circulation  at  the  trailing  edge  of  the  wing  section 
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The  boundary  condition  in  the  wake  is  to  impose  the  jump  in  potential, 

A(J>,  across  the  sheet.  When  the  wake  is  a  fixed,  planar  sheet,  the  jump  in  po¬ 
tential  is  held  constant  to  the  downstream  boundary  at  each  spanwise  station. 
For  a  free  wake  calculation,  the  jump  in  potential  is  constant  across  the  vor¬ 
tex  lines,  which  can  move  spanwise.  This  results  in  a  distribution  of  the 
jump  in  potential  that  varies  both  spanwise  and  streamwise. 

The  use  of  the  shearing  transformation  (Equation  11)  for  the  embedded 
fine  grid  system  complicates  the  imposition  of  symmetry  conditions  and  root 
juncture  conditions.  The  simple  Cartesian  (crude)  grid  symmetry  condition 
(Equation  3a)  becomes 


=  +  VY  =  ° 

at  the  symmetry  plane  and 

4 >Y  -  4>|Cy  +  4>nnY  -  Fx 


(25) 


(26) 


for  fuselage  combinations  in  the  root  juncture  region.  Here  represents  the 
slope  of  the  fuselage  or  body  at  the  root.  The  Z  shearing  is  not  used  because 
the  grid  is  always  continued  inside  the  root  as  a  constant  Z  value. 

Computations  indicate  that  numerical  instabilities  will  result  if  special 
attention  is  not  given  to  the  selection  of  difference  approximations  in  this 
region.  These  difficulties  result  from  the  nature  of  the  shearing  transforma¬ 
tion.  To  solve  this  problem,  a  plane  of  dummy  mesh  points  is  positioned 
across  the  symmetry  plane  or  within  the  body  surface.  These  flow-field  poten¬ 
tials  are  artificial  in  the  sense  that  there  is  no  physical  flow  field  associ¬ 
ated  with  them.  They  simply  provide  a  side  boundary  of  potentials  which,  when 
used  for  differencing,  produce  the  proper  side  condition  given  by  Equation 
(25)  or  (26). 

A  special  first  order  accurate  one-sided  difference  operator  is  used  to 
generate  the  dummy  interior  point  potential  values.  For  grid  lines  that  are 
swept  back  in  the  physical  plan  (^Y  <  0),  the  following  equation  is  used 
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and  for  grid  lines  that  are  swept  forward  (4^  >  0), 


D  /  []y  ^ Y 
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(27b) 


Note  that  the  operator  changes  depending  on  whether  the  grid  lines  are  swept 
forward  or  backward.  In  each  case,  th«»  coefficient  of  the  dummy  potential  at 
the  point  (I,  J-l)  is  larger  than  the  coefficients  of  other  potentials  in  the 
difference  equation.  This  enhances  the  effective  diagonal  dominance  of  the 
system  even  though  the  dummy  points  are  not  directly  relaxed  in  the  conven¬ 
tional  sense. 

Differencing  at  the  .  ing  tip  is  complicated  by  the  fact  that  the  fine 
mesh  system  does  not  extend  beyond  the  wing  tip.  Unlike  the  conventional  glo¬ 
bal  transformation  approach,  the  coordinate  lines  do  not  have  to  be  unswept  or 
unsheared  far  from  the  wing.  Provisions  must  be  made,  however,  for  properly 
ending  the  fine  grid  computation  at  the  wing  and  canard  tips.  For  this  rea¬ 
son,  another  temporary  fine  mesh  is  positioned  just  beyond  the  wing  tips. 

Like  its  neighboring  fine  grid  nearest  the  tip,  this  grid  array  is  located  at 
a  crude  mesh  Y-line.  Its  extent  in  the  streamwise  and  vertical  directions  is 
consistent  with  the  fine  grid  system.  Both  the  dummy  planes  beyond  the  root 
and  the  temporary  fine  grids  beyond  the  tips  are  computed  for  each  sweep  of 
thj  array  of  fine  grid  structures.  While  the  root  dummy  plane  is  computed 
using  difference  formulas,  the  tip  plane  is  simply  filled  using  linear  inter¬ 
polation  and  potentials  from  the  crude  Cartesian  grid. 


2.6  FUSELAGE  MODELING 

The  body  is  modeled  in  the  solution  by  a  constant  cross-section  computa¬ 
tional  surface  in  the  Cartesian  crude  grid.  The  input  data  allow  for  simple 
ax i symmetric  body  definition  or  detailed  "Quick  Geometry"  (Ref  3)  body  def¬ 
inition.  Body  boundary  conditions  are  imposed  by  fixing  the  velocity  poten- 
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Vachris,  A.  F.  and  Yaeger,  L.  S.,  "Quick-Geometry  -  A  Rapid  Response 
Method  for  Mathematically  Modeling  Configuration  Geometry,"  NASA  SP-390, 


tial  values  on  the  body  computational  surface.  The  procedure  follows  that  of 
Ref  2  and  is  described  in  detail  in  Ref  1.  Body  pressures  from  Equation  (2) 
are  used  to  produce  a  calculation  of  the  body  force  and  moment  contribution. 


2.7  BOUNDARY  LAYER  CALCULATION 

Viscous  effects  are  computed  in  the  analysis  code  by  coupling  a  modified 
Bradshaw  boundary  layer  computation  with  the  inviscid  potential  flow  solution. 
The  boundary  layer  calculation  is  virtually  identical  to  the  method  developed 
by  Mason  (Ref  4).  The  method  employs  the  modified  chord  technique  of  Nash 
(Ref  5),  which  represents  an  infinite  sheared  wing  boundary  layer  calculation. 
The  wing  sweep  angle  is  that  of  the  local  mid-chord  span  line,  such  that  it 
may  vary  across  the  span.  The  two  dimensional  Bradshaw  turbulent  boundary 
layer  analysis  (Ref  6)  provides  the  foundation  for  the  method.  The  use  of  a 
modified  two  dimensional  boundary  layer  analysis  greatly  reduces  the  necessary 
computer  time  and  has  demonstrated  good  results  for  several  different  codes 
(Ref  1,  2,  3,  7). 

The  boundary  layer  calculation  provides  a  displacement  thickness  and  skin 
friction  calculation  at  each  analysis  station  of  the  wing  and  canard.  The 
slope  of  the  displacement  thickness  is  used  to  modify  the  surface  boundary 
conditions  in  the  inviscid  solution.  The  local  skin  friction  calculation  is 
used  to  provide  a  viscous  drag  estimate  for  the  configuration  at  the  end  of 
the  analysis  run. 


4  Mason,  W.  H.,  et  al.,  "An  Automated  Procedure  for  Computing  the  Three- 
Dimensional  Transonic  Flow  Over  Wing-Body  Combinations,  including  Viscous 
Effects,"  Report  AFFDL-TR-77- 122,  Vol.  1,  October  1977. 

5  Nash,  J.  F.  and  Tseng,  R.  R.,  "The  Three-Dimensional  Turbulent  Boundary 
Layer  on  an  Infinite  Yawed  Wing,"  The  Aeronautical  Quarterly,  November 
1971. 

6  Bradshaw,  P.  and  Ferriss,  D.H.,  "Calculation  of  Boundary  Layer  Develop¬ 
ment  Using  the  Turbulent  Energy  Equation.  Compressible  Flow  on  Adiabatic 
Walls,"  J.  Fluid  Mech. ,  Vol.  46,  1971. 
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3  -  SAMPLE  CASES 


Two  fighter  configurations  and  one  transport  configuration  were  selected 
as  sample  cases  for  demonstration  of  the  CANTATA  code  multiple  surface  analy¬ 
sis  capability.  The  two  fighter  configurations  were  the  X-29  Forward  Swept 
Wing  technology  demonstrator  (Figure  4),  and  the  Configuration  Development  of 
Advanced  Fighters  (CDAF)  Configuration  (Figure  5).  The  transport  configura¬ 
tion  was  the  C-5A,  Figure  6. 

The  X-29  configuration  was  analyzed  at  Mach  0.9,  7.6  degrees  angle  of 
attack.  The  predicted  pressures  are  compared  with  wind  tunnel  data  (Ref  8)  in 
Figure  7.  Both  fixed  and  free  wake  analysis  results  are  shown.  The  discrep¬ 
ancies  between  the  analysis  and  data  are,  in  part,  due  to  flow  separation. 

The  test  Reynolds  number  was  approximately  four  million,  based  on  the  mean 
aerodynamic  chord.  The  computational  model  does  not  include  the  rearward 
strake  of  the  wind  tunnel  model,  and  the  code  does  not  allow  any  modeling  of 
inlet  spillage.  Lower  surface  actuator  fairings  were  included  in  the  wind 
tunnel  test,  but  not  modeled  in  the  analysis.  In  general,  the  code  predicts  a 
shock  location  that  is  aft  of  the  experimental  data,  except  at  0.907  semispan. 
There  the  agreement  with  data  is  very  good.  Considering  the  presence  of  the 
actuator  fairings  in  the  data,  the  lower  surface  agreement  is  very  good.  The 
upper  surface  leading  edge  expansion  is  under-predicted  outboard  of  0.306 
semispan.  The  results  at  0.907  semispan  might  be  considered  typical  of  a 
small  disturbance  prediction,  but  the  discrepancies  at  0.490  and  0.698  are 
much  greater.  Model  inspection  results  were  not  available  to  verify  the 
leading  edge  contours.  The  free  wake  analysis  produces  a  more  forward  shock 

8  Charletta,  Roy,  "Post  Test  Report  Series  I  Transonic/Supersonic  Testing 

on  a  12.5’A  Scale  Grumman  Design  712,  X-29A  Forward  Swept  Wing  Demonstrator 
Aircraft  Model  in  the  NASA-ARC  11  foot  and  9x7  foot  Wind  Tunnels,  at 
Moffett  Field,  CA,"  Grumman  Aerospace  Report  No.  712/ENG-RPT82-021 , 

August  1982. 
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Figure  6.  C-SA  General  arrangement. 
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Figure  7.  Date-analysis  pressure  comparison,  X-29,  Mach  0.9,  a  7.6  (sheet  2  of  2). 


location  inboard,  with  very  little  effect  outboard.  This  would  be  expected, 
considering  that  the  wing  tip  is  fairly  well  separated,  both  outboard  and  up¬ 
stream,  from  the  principal  canard  wake  effects. 

Analysis  results  for  the  CDAF  configuration  at  Mach  0.9,  7.0  degrees 
angle-of-attack  are  compared  with  data  (Ref  9)  in  Figure  8.  The  flow  condi¬ 
tions  are  severe  for  a  potential  flow  analysis,  with  local  Mach  numbers  ex¬ 
ceeding  1.7  at  several  locations.  The  canard  and  inboard  wing  pressures  show 
good  agreement  with  the  data.  The  two  outboard  wing  stations  are  not  as  well 
predicted.  This  may  be  because  the  local  flow  conditions  are  too  severe  for 
potential  flow  modeling,  or  because  of  flow  separation  which  is  beyond  the 
modeling  capability  of  the  boundary  layer  analysis.  The  free  wake  analysis 
results  in  a  more  negative  leading  edge  pressure  peak  at  the  two  inboard  sta¬ 
tions.  Relative  to  a  wake  shape  that  is  fixed  at  the  height  of  the  canard 
surface,  a  "free"  wake  moves  downward  (for  positive  canard  lift).  The  results 
shown  in  Figure  8  indicate  that  a  more  complex  process  occurs.  Figure  9  shows 
the  spanwise  movement  of  the  canard  and  wing  vortex  lines.  There  is  a  sub¬ 
stantial  redistribution  of  the  vorticity  in  the  canard  wake.  The  edge  of  the 
canard  wake  moves  inboard  more  than  one  grid  plane,  and  two  inboard  canard 
vortex  lines  merge  together  just  ahead  of  streamwise  station  100. 

Force  and  moment  predictions  for  the  CDAF  conf iguration  are  presented  in 
Figure  10.  Three  angles-of-attack,  6.0,  7.0  and  8.0  degrees,  and  three  canard 
settings,  0,  -5  and  -10  degrees,  were  analyzed.  The  comparison  with  data 
shows  good  prediction  of  stability  level  and  canard  effectiveness.  An  overall 
shift  in  pitching  moment  of  +0.06  has  been  added  to  the  predictions.  Several 
effects  acting  together  are  the  most  likely  cause  of  this  shift.  The  under¬ 
prediction  of  the  leading  edge  expansion  would  result  in  a  more  negative 
pitching  moment  prediction.  In  addition,  the  prediction  of  more  negative  pres¬ 
sures  near  the  upper  surface  trailing  edge  and  more  positive  pressures  in  the 
lower  surface  cove  result  in  more  negative  pitching  moment.  Perhaps  most  sig¬ 
nificant  is  the  lack  of  a  detailed  body  solution.  The  body  modeling  is  ade¬ 
quate  to  account  for  induced  effects  on  the  lifting  surfaces,  but  is  known  to 

9  Spurlin,  C.  J.,  "Test  Report  for  AFWAL  Optimal  Transonic  Configuration 

Fighter,  Project  No.  P41T-G6,  Test  No.  TF-570,"  April  1980. 
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Figure  10.  Fores  and  moment  prediction,  CDAF,  Mach  0.9 
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be  inadequate  near  the  fuselage  nose.  Thus,  the  predicted  body  contribution, 
to  pitching  moment  is  unreliable.  The  free  wake  modeling  results  in  a 
slightly  more  negative  pitching  moment  and  a  more  stable  pitching  moment 
slope. 

The  predicted  drag  is  compared  to  wind  tunnel  data  in  Figure  10b.  The 
predicted  drag  levels  for  the  fixed  wake  analysis  are  all  too  high.  As  for 
the  pitching  moment  prediction,  the  lack  of  an  accurate  body  solution  is  prob¬ 
ably  part  of  the  error.  In  general,  a  small  disturbance  analysis  is  not  ex¬ 
pected  to  predict  the  leading  edge  suction  of  the  lifting  surfaces.  Thus,  the 
increase  in  drag  for  an  increase  in  lift  will  be  overpredicted,  as  the  results 
show.  The  free  wake  analysis  results  in  lower  predicted  drag  at  all  lift  co¬ 
efficients.  At  -  0.5,  the  predicted  drag  values  are  very  good,  falling 
within  the  experimental  uncertainty  (Ref  1).  As  for  the  fixed  wake  analysis, 
the  increase  in  drag  for  an  increase  in  lift  is  overpredicted. 

The  C~5A  represents  a  much  milder  case  for  potential  flow  analysis.  Fig¬ 
ure  11  compares  data  (Ref  10)  and  analysis  results  at  Mach  0.75  and  2.0 
degrees  angle-of-attack.  In  general,  the  leading  edge  expansion  on  the  wing 
is  under-predicted,  and  the  predicted  outboard  shock  is  too  far  aft.  The 
predicted  horizontal  tail  pressures  show  good  agreement  with  the  data.  The 
inboard  station  on  the  horizontal  tail  is  influenced  by  the  "bullet"  at  the 
top  of  the  vertical  tail.  The  bullet  is  not  modeled  in  the  analysis.  The 
effect  of  the  free  wake  model  is  slight  on  both  the  wing  and  tail,  slightly 
improving  the  agreement  on  the  horizontal  tail. 

During  the  study,  pressure  data  for  a  "Propulsive  Wing-Canard"  configura¬ 
tion  became  available  (Ref  11).  The  configuration  is  shown  in  Figure  12,  and 
the  comparison  of  predicted  and  experimental  pressures  is  shown  in  Figure  13. 
Generally  good  agreement  is  seen,  with  the  data  indicating  trailing  edge  sepa¬ 
ration  much  more  than  the  analysis. 
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10  Harris,  M.  K. ,  Huie,  W.  E,  "C-5A  Aerodynamic  Data  for  Airloads," 
Lockheed-Georgia  Report  LG74ER0162,  November  1974. 

11  Stewart,  V.  R.,  "Evaluation  of  a  Propulsive  Wing/Canard  Concept  at 
Subsonic  and  Supersonic  Speeds,"  Naval  Air  Systems  Command  Report 
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d)  Wing,  2y/b  =  0.705 


Figure  13.  Data-analysis  pressure  comperison,  propulsive  wing-canard,  Mach  0.9,  a  4.0  (sheet  2  of  2) 


4  -  COPES/CONMIN  OPTIMIZATION 


The  PANDORA  wing-body-canard  analysis  code  was  coupled  with  the  COPES  and 
CONMIN  routines  of  Vanderplaats  (Ref  12,  13).  The  COPES  code  is  a  control 
program  that  connects  the  numerical  optimization  code  CONMIN  with  the  aero¬ 
dynamic  analysis  code.  The  parts  of  the  PANDORA  code  allowing  the  numerical 
optimization  have  been  left  intact  in  the  CANTATA  code. 

The  COPES  and  CONMIN  routines  were  slightly  modified  in  the  PANDORA  code. 
Changes  were  made  in  the  main  COPES  routine  and  subroutine  CNM06 .  The  change 
to  CNM06  is  important.  It  provides  subroutine  ANALIZ  with  the  index  of  the 
best  result  during  an  optimization  search.  In  this  way,  the  next  search  solu¬ 
tion  can  be  restarted  from  the  best  previous  result.  This  more  closely  models 
the  search  strategy  in  CONMIN,  resulting  in  more  consistent  information  from 
the  aerodynamic  analysis.  The  changes  in  the  COPES  and  CONMIN  routines  are 
identified  in  the  code  with  comment  cards. 

The  optimization  algorithm  within  CONMIN  is  a  modified  Method  of  Feasible 
Directions  (Ref  13).  The  gradient  information  for  the  algorithm  is  calculated 
by  sequentially  perturbing  each  design  variable.  Each  design  variable  per¬ 
turbation  requires  analysis  by  the  flow  solution  routines.  Thus,  the  computer 
time  is  proportional  to  the  number  of  design  variables.  The  gradient  informa¬ 
tion  establishes  a  search  direction  that  should  improve  the  design  (decrease 
the  objective  function)  while  satisfying  any  constraints.  The  search  direc¬ 
tion  is  "explored"  until  a  relative  optimum  is  found  or  any  constraints  become 
violated.  If  the  starting  conditions  violated  any  constraints,  the  search 
direction  will  be  that  which  satifies  the  constraints  with  the  least  objective 
function  increase.  One  to  four  flow  solutions  are  required  during  the  search. 

12  Vanderplaats,  G.  N. ,  "COPES — Control  Program  for  Engineering  Synthesis," 
to  be  published  as  a  Naval  Postgraduate  School  Memorandum. 

13  Vanderplaats,  G.  N.,  "CONMIN — A  Fortran  Program  for  Constrained  Function 
Minimization,"  NASA  TM  X-62282,  August  1973. 


Completion  of  the  search  constitutes  one  optimization  iteration.  The  follow¬ 
ing  discussion  of  COPES  usage  is  excerpted  from  Ref  12.  A  more  detailed  de¬ 
scription  will  be  found  in  the  reference.  Reference  13  should  be  consulted  if 
more  details  of  the  numerical  optimization  algorithm  are  desired. 

If  it  is  desired  to  run  only  a  simple  analysis  using  COPES,  only  three 
data  cards  are  required  for  the  COPES  program:  a  TITLE  card,  a  control 
parameter  (NCALC  =  1),  and  an  END  card.  If  the  optimization  or  parametric 
analysis  (sensitivity)  capabilities  of  COPES  are  to  be  used,  additional  data 
must  be  read.  This  data  will  identify  which  parameters  in  the  global  common 
block,  GLOBCM,  are  used.  To  set  up  the  COPES  data,  the  user  must  have  a  basic 
understanding  of  how  the  data  in  the  global  common  block  is  accessed  by  COPES. 
This  is  outlined  in  the  following  section. 

4.1  OPTIMIZATION  DATA  MANAGEMENT 

In  order  to  perform  design  operations,  the  COPES  program  must  access  the 
data  in  common  block  GLOBCM.  This  is  done  by  defining  the  location  in  GLOBCM 
where  a  specified  parameter  resides.  For  example,  consider  the  common  block 
for  a  cantilevered  beam  design  problem: 

COMMON/GLOBCM/B , H , VOL , BSTRES , SHRSTR , DELTA , HB , E , AL 

The  volume  of  material,  VOL,  is  the  third  parameter  in  the  common  block; 
that  is,  it  resides  in  Location  3,  referred  to  as  the  Global  Location  number. 
Similarly  the  bending  stress,  BSTRES,  is  in  Global  Location  4  and  the  beam 
width  is  in  Global  Location  1.  Thus,  the  parameters  are  referred  to  by  their 
respective  location  numbers  in  global  common. 

For  convenience  in  preparing  data  for  the  COPES  program,  a  simple 
"CATALOG"  of  parameters  may  be  defined.  For  the  cantilevered  beam,  this 
catalog  would  be: 


GLOBAL 

FORTRAN 

LOCATION 

NAME 

DEFINITION 

1 

B 

Beam  width 

2 

H 

Beam  height 

3 

VOL 

Volume  of  material 

4 

BSTRES 

Maximum  bending  stress 

5 

SHRSTR 

Maximum  shear  stress 

6 

DELTA 

Deflection  under  the  load 

7 

HB 

Ratio,  H/B 

8 

E 

Young's  modulus 

9 

AL 

Length  of  beam 

As  another  example,  consider 

a  global  common 

block  containing  arrays: 

GLOBAL 

FORTRAN 

LOCATION 

NAME 

DEFINITION 

1 

A 

Area 

2 

Y(  10) 

Vector 

12 

Q 

a 

13 

C(2,2) 

• 

17 

H 

etc. 

The  dimensions  are  given  with  the  FORTRAN  name  as  a  reminder  that  the 
parameter  is  an  array.  In  this  case,  the  third  parameter  in  the  Y  array  is  in 
Global  Location  4.  Remembering  that  arrays  are  stored  column  by  column,  the 
C(l,2)  array  location  is  in  Global  Location  15. 

It  will  be  seen  that  identifying  parameters  according  to  their  location 
in  GLOBCM  provides  a  great  deal  of  flexibility  in  using  the  COPES  program  for 
des ign. 

In  the  following  section,  definitions  of  terms  commonly  used  in  automated 
design  are  given  for  easy  reference. 

4.2  OPTIMIZATION  TERMINOLOGY 


The  COPES  program  currently  provides  six  specific  capabilities: 
•  Simple  analysis  -  just  as  if  COPES  was  not  used 


•  Optimization  -  minimization  or  maximization  of  one  calculated  function 
Kith  limits  imposed  on  other  functions 

•  Sensitivity  analysis  -  the  effect  of  changing  one  or  more  design  vari¬ 
ables  on  one  or  more  calculated  functions 

e  Two-variable  functions  space  -  analysis  for  all  specified  combinations 
of  two  design  variables 

e  Optimum  sensitivity  -  same  as  sensitivity  analysis  except,  at  each 
step,  the  design  is  optimized  with  respect  to  the  remaining  independ¬ 
ent  design  variables 

•  Approximate  optimization  -  optimization  using  approximation  tech¬ 
niques.  Usually  more  efficient  than  standard  optimization  for  up  to 
10  design  variables  if  multiple  optimizations  are  to  be  performed. 

In  defining  the  data  required  to  execute  the  COPES  program,  the  following 
definitions  are  useful. 

DESIGN  VARIABLES  -  Those  parameters  which  the  optimization  program  is 
allowed  to  change  in  order  to  improve  the  design.  Design  variables  appear 
only  on  the  right  hand  side  of  equations  in  the  analysis  program.  COPES 
considers  two  types  of  design  variables,  independent  and  dependent.  If  two  or 
more  variables  are  always  required  to  have  the  same  value  or  be  in  a  constant 
ratio,  one  is  the  independent  variable  while  the  remaining  are  dependent 
variables.  For  example,  if  the  height  is  required  to  be  10  times  the  width  of 
the  cantilevered  beam,  B  would  be  the  independent  variable  while  H  would  be 
the  dependent  variable. 

OBJECTIVE  FUNCTION  -  The  parameter  to  be  minimized  or  maximized  during 
optimization.  Also  the  parameters  calculated  as  functions  of  specified  design 
variables  during  a  sensitivity  or  two-variable  function  space  study. 

Objective  functions  always  occur  on  the  left  side  of  equations,  unless  the 
objective  function  is  also  a  design  variable  (the  beam  height  may  be  minimized 
as  an  objective  function  if  it  is  also  a  design  variable;  in  this  way,  the 
minimum  height  is  found  for  which  no  constraints  are  violated).  An  objective 
function  may  be  linear  or  non-linear,  implicit  or  explicit,  but  must  be  a 
function  of  the  design  variables  to  be  meaningful. 

CONSTRAINT  -  Any  parameter  which  must  not  exceed  specified  bounds  for  the 
design  to  be  acceptable.  Constraint  functions  always  appear  on  the  left  side 
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of  equations.  Just  as  for  objective  functions,  constraints  may  be  linear  or 
non-linear,  implicit  or  explicit,  but  must  be  functions  of  the  design  variables, 

CONSTRAINT  SET  -  A  group  of  constraints  which  appear  consecutively  in  the 
global  common  block  and  which  all  have  the  same  limits  imposed.  This  is  a 
convenience  which  allows  several  constraints  to  be  identified  with  a  minimum 
of  data. 

GLOBAL  COMMON  -  Common  block  GLOBCM  containing  design  information. 

GLOBAL  LOCATION  -  Location  of  a  particular  parameter  in  GLOBCM. 

4.3  OPTIMIZATION  INPUT  DATA  FORMAT 

In  order  to  execute  the  COPES  program  it  is  necessary  to  provide 
formatted  data  for  COPES,  followed  by  data  for  the  ANALIZ  program  which  is 
coupled  to  COPES.  Section  4.1  defines  the  data  which  is  required  by  COPES. 

The  data  is  segmented  into  "BLOCKS"  for  convenience.  All  formats  are  alpha¬ 
numeric  for  TITLE,  and  END  cards,  FIO  for  real  data  and  110  for  integer  data. 

The  COPES  data  begins  with  a  TITLE  card  and  ends  with  an  END  card.  This 
is  followed  by  data  to  be  read  by  the  user-supplied  subroutine  ANALIZ. 

Comment  cards  may  be  inserted  anywhere  in  the  COPES  data  stack  prior  to 
the  END  card,  and  are  identified  by  a  dollar  sign  ($)  in  column  1. 

While  the  input  description  defines  COPES  data  in  formatted  fields  of 
ten,  the  data  may  actually  be  read  in  more  conveniently  by  separating  data  by 
commas  or  one  or  more  blanks.  If  more  than  one  number  is  contained  on  an  un¬ 
formatted  data  card,  a  comma  must  appear  somewhere  on  the  card.  If  exponen¬ 
tial  numbers  such  as  2.5+10  are  read  on  an  unformatted  card,  there  must  be  no 
embedded  blanks.  Unformatted  cards  may  be  intermingled  with  formatted  cards. 
Real  numbers  on  an  unformatted  card  must  have  a  decimal  point.  Some  examples: 

Unformatted  data: 

5,7, 1.3, 1.0+20,0, -5.1 
5,7,1.3,1.0+20, ,-5.1 
5  7  1.3  1.0+20,,  -  5.1 
5  7  1.3,  1.0+20  0  -5.1 

Equivalent  formatted  data: 

col _ 10 _ 20 _ 30 _ 40 _ 50 _ 60 _ 70 _ 80 

5  7  1.3  1.0+20  0  -5.1 


Unformatted  data 
2 

2.3 
2  3 

Equivalent  formatted  data: 

col _ 10  20  30 _ 

2 

2  3 

2  3 

Note  that  this  data  contains  no  commas,  so  it  is  assumed  to  be  formatted  al¬ 
ready. 

Unformatted  data: 

1.2.3.4.5.6.7.8.9.10.11 

col _ 10 _ 20 _ 30  40 _ 50 _ 60 _ 70 _ 80 

1  2  3  4  5  6  7  8 

9  10  11 

Note  that  two  formatted  data  cards  are  created  here. 

Unformatted  data: 

1,2, 3, 4, 5, 6 

7.8.9.10.11 

Equivalent  formatted  data: 

col _ 10 _ 20  30 _ 40 _ 50 _ 60 _ 70 _ 80 

1  2  3  4  5  6 

7  8  9  1C  11 

Note  that  the  two  examples  above  do  not  produce  the  same  formatted  data  cards. 
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5  -  INPUT  DATA  DESCRIPTION 


This  section  describes  the  necessary  solution/optimization  parameters  and 
geometry  input  needed  to  operate  the  CANTATA  code.  The  information  of  Section 
4,  COPES/CONMIN  Optimization,  is  important  for  using  the  optimization.  Only 
the  one-cycle  analysis  (NCALC  *  1  in  DATA  BLOCK  B)  and  optimization  (NCALC  « 

2)  COPES  options  have  been  thoroughly  tested  in  this  effort.  The  transonic 
analysis  input  is  an  extension  of  that  from  the  base  code  (Ref  1).  A  set  of 
Usage  Notes  (Subsection  5.3)  appears  at  the  end  of  the  analysis  input 
description  to  elaborate  on  the  less  familiar  input  parameters.  Volume  2, 

Part  2  of  Ref  1  contains  important  discussions  of  the  optimization  application 
that  are  not  repeated  here. 

5.1  OPTIMIZATION  INPUT  DATA 

DATA  BLOCK  A_ 

DESCRIPTION:  Title  card. 

FORMAT  AND  EXAMPLE 

12  3  4  5  6  7  8 

TITLE 

CANTILEVERED  BEAM  DESIGN 

FIELD  CONTENTS 

1-8  Any  80  character  title  may  be  given  on  this  card. 


FORMAT 

20A4 


DATA  BLOCK 
DESCRIPTION: 

FORMAT  AND  EXAMPLE 


B 

Program  Control  Parameters. 


1 

2  3 

4 

5 

6 

7 

FORMAT 

NCALC 

NDV  NSV 

N2VAR 

NXAPRX 

IPNPUT 

IPDBG 

7110 

2 

2  3 

5 

2 

0 

0 

FIELD 

CONTENTS 

1  NCALC 

:  Calculation  Control 

0  -  Read  input  and  stop.  Data  of  blocks  A,  B  and  V  is  required. 

Remaining  data  is  optional. 

1  -  One  cycle  through  program.  The  same  as  executing  ANALIZ 

stand-along  (i.e.,  no  optimization).  Data  of  blocks  A,  B  and 
V  is  required.  Remaining  data  is  optional. 

2  -  Optimization.  Data  of  blocks  A-I  and  V  is  required. 

Remaining  data  is  optional. 

3  -  Sensitivity  analysis.  Data  of  blocks  A,  B,  P,  Q  and  V  is 

required.  Remaining  data  is  optional. 

4  -  Two  variable  function  space.  Data  of  blocks  A,  B,  and  R-V  is 

required.  Remaining  data  is  optional. 

5  -  Optimum  sensitivity.  Data  of  blocks  A-I,  P,  Q,  and  V  is 

required.  Remaining  data  is  optional. 

6  -  Optimization  using  approximation  techniques.  Data  of  blocks 

A-0  and  V  is  required.  Remaining  data  is  optional. 

2  NDV :  Number  of  design  variables  on  which  sensitivity  analysis  or 

optimization  will  be  performed. 

4  N2VAR :  Number  of  objective  functions  in  a  two  variable  function 

space  study. 

5  NXAPRX :  Number  of  X-variables  for  approximate  analysis/optimization. 

>  NDV 

6  IPNPUT:  Input  print  control 

0  -  Print  card  images  of  data  plus  formatted  print  of  input  data 

1  -  Formatted  print  only  of  input  data 

2  -  No  print  of  input  data. 


DATA  BLOCK  OMIT  IF  NDV  =  0  IN  BLOCK  B 

DESCRIPTION:  Integer  optimization  control  parameter. 


FORMAT  AND 

1 

EXAMPLE 

2 

3 

4 

5 

6 

7 

8 

FORMAT 

IPRINT 

ITMAX 

ICNDIR 

NSCAL 

ITRM 

LINOBJ 

NACMXl 

NFDG 

7110 

5 

0 

0 

0 

0 

0 

0 

0 

FIELD 

CONTENTS 

1  IPRINT:  Print  control  used  in  the  optimization  program  CONMIN. 

0  -  No  print  during  optimization. 

1  -  Print  initial  and  final  optimization  information. 

2  -  Print  above  plus  objective  function  value  and  design  variable 

values  at  each  iteration. 

3  -  Print  above  plus  constraint  values,  direction  vector  and  move 

parameter  at  each  iteration. 

4  -  Print  above  plus  gradient  information. 

5  -  Print  above  plus  each  proposed  design  vector,  objective 

function  and  constraint  values  during  the  one-dimensional 
search. 

2  UMAX:  Maximum  number  of  optimization  iterations  allowed. 

DEFAULT  -  20. 

3  ICNDIR :  Conjugate  direction  restart  parameter.  DEFAULT  -  NDV  +  1. 

4  NSCAL :  Scaling  parameter.  .GT.O  -  Scale  design  variables  to  order 

of  magnitudes  one  every  NSCAL  iterations.  .LT.O  -  Scale 
design  variables  according  to  user-input  scaling  values. 

Good  value  is  ICNDIR,  1  is  not. 

5  ITRM :  Number  of  consecutive  iterations  which  must  satisfy  relative 

or  absolute  convergence  criterion  before  optimization  process 
is  terminated.  DEFAULT  -  3 

6  LINOBJ:  Linear  objective  function  identifier.  If  the  optimization 

objective  is  known  to  be  a  linear  function  of  the  design 
variables,  set  LINOBJ  “  1.  DEFAULT  -  Nonlinear. 

7  NACMXl:  One  plus  the  maximum  number  of  active  constraints 

anticipated.  DEFAULT  -  NDV  +  2. 


8 


NFDG : 


Finite  difference  gradient  identifier. 

0  -  All  gradient  information  is  computed  by  finite  difference 

with  CONMIN. 

1  -  All  gradient  information  is  computed  analytically  by  the 

user-supplier  code. 

2  -  Gradient  of  objective  is  computed  analytically.  Gradients  of 

constraints  are  computed  by  finite  difference  within  CONMIN.' 

REMARKS 

•  Currently  NFDG  must  be  zero  in  COPES 

•  IPRINT  =  5  is  recommended. 


i 

i 

i 

I 


( 

i 


j 


I 
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DATA  BLOCK 

£  OMIT 

IF  NDV  »  0  IN  BLOCK 

B 

DESCRIPTION 

:  Floating 

point  optimization 

program  parameters 

• 

FORMAT  AND 

EXAMPLE 

1 

2 

3  4 

5 

6 

7 

FORMAT 

FDCH 

FDCHM 

CT  CTMIN 

CTL 

CTLMIN 

THETA 

7F10 

0.0 

0.0 

0.0  0.0 

0.0 

0.0 

0.0 

DELFUN 

DAB FUN 

ALPHAX  ABOBJ1 

0.0 

0.0 

0.0  0.0 

NOTE: 

Two  cards 

are  read  here. 

FIELD 

CONTENTS 

1 

FDCH: 

Relative  change  in  design  variables  in  calculating  finite 

difference  gradients.  DEFAULT  *  0.01. 

2 

FDCHM : 

Minimum  absolute  step  in  finite  difference  gradient 

calculations.  DEFAULT  =■  0.001. 

3 

CT: 

Constraint  thickness  parameter.  DEFAULT  -  0.05. 

4 

CTMIN: 

Minimum  absolute  value  of  CT  considered  in  the  optimization 

process.  DEFAULT  =  0.004. 

5 

CTL: 

Constraint  thickness  parameter  for  linear  constraints. 

DEFAULT  =  -0.01. 

6 

CTLMIN: 

Minimum  absolute  value  of  CTL  considered  in  the  optimization 

process.  DEFAULT  *  0.001. 

7 

THETA: 

Mean  value  of  push-off  factor  in  the  method  of  feasible 

directions.  DEFAULT  =  1.0. 

1 

DELFUN: 

Minimum  relative  change  in  objective  function  to  indicate 

convergence  of  the  optimization  process.  DEFAULT  •  0.001. 

2 

DABFUN : 

Minimum  absolute  change  in  objective  function  to  indicate 

convergence  of  the  optimization  process.  DEFAULT  -  0.001. 

3 

ALPHAX: 

Maximum  fractional  change  in  any  design  variable  for  first 

estimate  of  the  step  in  the  one-dimensional  search. 

DEFAULT  -  0.1. 

4 

ABOBJ1 : 

Expected  fractional  change  in  the  objective  function  for 

first  estimate  of  the  step  in  the  one-diroensional  search. 
DEFAULT  -  0.1. 


REMARKS 


•  The  DEFAULT  values  for  these  parameters  usually  work  well. 


DATA  BLOCK 


E  OMIT  IF  NDV  -  0  IN  BLOCK  B 


DESCRIPTION:  Total  number  of 

design  variables,  design  objective 

identification 

and  sign. 

FORMAT  AND  EXAMPLE 

1  2 

3 

FORMAT 

NDVTOT  IOBJ 

SGNOPT 

2110, F10 

0  3 

-1.0 

FIELD 

CONTENTS 

1  NDVTOT :  Total  number  of  variables  linked  to  the  design  variables. 

This  option  allows  two  or  more  parameters  to  be  assigned  to  a 
single  design  variable.  The  value  of  each  parameter  is  the 
value  of  the  design  variable  times  a  multiplier,  which  may  be 
different  for  each  parameter.  DEFAULT  *  NDV. 

2  IOBJ:  Global  variable  location  associated  with  the  objective  func¬ 

tion  in  optimization. 

3  SGNOPT:  Sign  used  to  identify  whether  function  is  to  be  maximized  or 

minimized.  +1.0  indicates  maximization.  -1.0  indicates  min¬ 
imization.  If  SGNOPT  is  not  unity  in  magnitude  it  scales  the 
magnitude  of  the  objective. 
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DATA  BLOCK  JF  OMIT  IF  NDV  *=  0  IN  BLOCK  B 

DESCRIPTION:  Design  variable  bounds,  initial  values  and  scaling  factors. 


FORMAT  AND  EXAMPLE 


FORMAT 


5:5 

A* 

A* 

,IT*' 

n 


.5  5.  0.0  0.0 

NOTE:  Read  one  card  for  each  of  the  NDV  independent  design  variables. 

Values  are  in  order  for  design  variables  NDSGN  sequence  (i.e.,  input 
sequence) . 


FIELD 


SCAL: 


CONTENTS 

Lower  bound  on  the  design  variable.  If  VLB.LT.-1.0E+15,  no 
lower  bound. 

Upper  bound  on  the  design  variable.  If  VUB.GT. 1.0E+15,  no 
upper  bound. 

Initial  value  of  the  design  variable.  IF  X  is  non-zero,  this 
will  supersede  the  value  initialized  by  the  user-supplied 
subroutine  ANALIZ. 

Design  variable  scale  factor.  Not  used  if  NSCAL.GE.O  in 
BLOCK  C. 


DATA  BLOCK  £  OMIT  IF  NDV  -  0  IN  BLOCK  B 
DESCRIPTION:  Design  variable  identification. 

FORMAT  AND  EXAMPLE 

1  2  3  FORMAT 

NDSGN  IDSGN  AMULT  2I10.F10 

1  1  1.0 


NOTE: 

FIELD 

1 

2 

3 


Read  one  card  for  each  of  the  NDVTOT  design  variables. 

CONTENTS 

NDSGN:  Design  variable  number  associated  with  this  variable. 

IDSGN:  Global  variable  number  associated  with  this  variable. 

AMULT:  Constant  multiplier  on  this  variable.  The  value  of  the 

variable  will  be  the  value  of  the  design  variable,  NDSGN, 
times  AMULT.  DEFAULT  -  1.0. 


DATA  BLOCK  H  OMIT  IF  NDV  -  0  IN  BLOCK  B 
DESCRIPTION:  Number  of  constrained  parameters. 

FORMAT  AND  EXAMPLE 

1  FORMAT 

NCONS  110 

4 

FIELD  CONTENTS 

1  NCONS:  Number  of  constraint  sets  in  the  optimization  problem. 

REMARKS 

•  If  two  or  more  adjacent  parameters  in  the  global  common  block  have  the 
same  limits  imposed,  these  are  part  of  the  same  constraint  set. 
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DATA  BLOCK  _I_  OMIT  IF  NDV  -  0  IN  BLOCK  B,  OR  NCONS  -  0  IN  BLOCK  H 
DESCRIPTION:  Constraint  identification  and  constraint  bounds. 

FORMAT  AND  EXAMPLE 


1  1 

2 

3 

4 

FORMAT 

|  ICON 

JCON 

LCON 

3110 

!  4 

0 

0 

BL 

SCAL1 

BU 

SCAL2 

-1.0+20 

0.0 

20000. 

0.0 

4F10 

i 


NOTE:  Read  two  cards  for  each  of  the  NCONS  constraint  sets. 

FIELD  CONTENTS 


1 

ICON: 

First  global  number  corresponding  to  the  constraint  set. 

2 

JCON: 

Last  global  number  corresponding  to  the  constraint  set. 

DEFAULT  -  ICON. 

3 

LCON: 

Linear  constraint  identifier  for  this  constraint  set.  LCON  ” 

1  indicates  linear  constraints.  If  in  doubt,  use  nonlinear. 

1 

BL: 

Lower  bound  on  the  constrained  variables.  If  BL.LT.-1.0E+15, 

no  lower  bound. 

2 

SCALl : 

Normalization  factor  on  lower  bound.  DEFAULT  ■  MAX  of 

(abs(bl) ,  0.1). 

3 

BU: 

Upper  bound  on  the  constrained  variables.  If  BU.GT. 1.0E+15, 

no  upper  bound. 

4 

SCAL2: 

Normalization  factor  on  upper  bound.  DEFAULT  -  MAX  of 

(ABS(BU),  0.1). 

REMARKS 

•  The  normalization  factor  should  usually  be  defaulted 

•  The  constraint  functions  sent  to  CONMIN  are  of  the  form:  (BL  - 
VALUE )/SCALl  .LE.  0.0  AND  (VALUE  -  BU)/SCAL2  .LE.  0.0 

•  Each  constrained  parameter  is  converted  to  two  constraints  in  CONMIN 
unless  ABS(BL)  or  ABS(BU)  exceeds  1.0E+15,  in  which  case  no  constraint 
is  created  for  that  bound. 


DATA  BLOCK  _J_  OMIT  IF  NXAPRX  -  0  IN  BLOCK  B 

DESCRIPTION:  Approximate  analysis/optimization  control  parameters. 
FORMAT  AND  EXAMPLE 


1 

2 

3 

A 

5 

6 

7 

8 

FORMAT 

NF 

NXS 

NXFS 

NXA 

INOM 

ISCRX 

IXCRXF 

IPAPRX 

8110 

5 

1 

1 

1 

0 

0 

0 

1 

KMIN 

KMAX 

MPMAX 

JNOM 

INXLOC 

INFLOC 

6110 

0 

0 

0 

0 

0 

0 

FIELD  CONTENTS 

1  NF:  Number  of  functions  to  be  approximated.  DEFAULT  -  number  of 

optimization  objective  and  constraint  functions. 

2  NXS:  Number  of  X-vectors  read  as  data. 

3  NXFS:  Number  of  X-F  pairs  read  as  data. 

A  NXA:  If  non-zero,  the  design  variables  read  by  SUBROUTINE  ANALIZ 

form  an  X-vector. 

5  INOM:  Nominal  X-vector  about  which  to  do  Taylor  expansion.  DEFAULT 

=  best  available. 

6  ISCRX:  File  from  which  NXS  X-vectors  are  read.  DEFAULT  •  5. 

7  IXCRXF:  File  from  which  NXFS  X-F  pairs  of  data  are  read.  DEFAULT  - 

5. 

8  IPAPRX:  Print  control,  1  to  A.  A  is  most. 

1  KMIN :  Minimum  number  of  approximation  iterations. 

2  KMAX:  Maximum  number  of  approximation  iterations. 

3  NPMAX:  Maximum  number  of  designs  retained  for  Taylor  series 

expansion. 

A  JNOM:  Number  of  iterations  after  which  the  best  design  is  picked  as 

nominal . 

5  INXLOC:  X-variable  global  location  identifier.  If  INXLOC  -  0,  the 

Taylor  series  expansion  is  on  the  design  variables  listed  in 
BLOCK  G. 

6  INFLOC:  Function  global  location  identifier.  If  INFLOC  *  0,  the 

Objective  and  constraint  functions  identified  in  BLOCKS  E  rnd 
I  are  the  functions  on  which  the  Taylor  series  expansion  is 
performed. 


REMARKS 


•  If  ISCRX  and/or  ISCRXF  file  number  is  other  than  5,  the  data  read  from 
that  file  is  assumed  to  be  binary  data 

•  If  NXS  -  NSFS  -  0,  NXA  is  defaulted  to  NXA  -  1,  even  if  it  is  read  as 
zero.  Also,  a  second  vector  of  design  variables  is  automatically 
defined  by  COPES  to  yield  two  independent  designs  to  start  the 
optimization. 
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DATA  BLOCK  K  OMIT  IF  NDV  -  0  IN  BLOCK  B,  OR  NXAPRX  -  0  IN  BLOCK  B 
DESCRIPTION:  Bounds  and  multipliers  for  approximate  optimization. 
FORMAT  AND  EXAMPLE 


DX2  DX3  DX4  DX5 


XFACT1  XFACT2 


8  FORMAT 


i 

'i1 


FIELD 


Two  or  more  cards  are  read  here. 

CONTENTS 

DXI:  Allowable  change  (in  magnitude)  of  the  Ith  design  variable 

during  each  approximate  optimization. 

XFACTl :  Multiplier  on  DXI  when  the  diagonal  elements  of  the  H  matrix 
are  available.  DEFAULT  -  1.5. 

XFACT2:  Multiplier  on  DFXI  when  all  elements  of  the  H  matrix  are 
available.  DEFAULT  =  2.0. 


DATA  BLOCK  L  OMIT  IF  NXAPRX  -  0  IN  BLOCK  B  OR  INXLOC  -  0  IN  BLOCK  J 
DESCRIPTION:  Global  locations  of  approximating  variables. 


FORMAT  AND  EXAMPLE 


LOCXl  LOCX2  LOCX3  LOCX4 


FORMAT 


•5 

•i 

'5 

$ 

•!i 


l/J 

I 

I L 

f 

i 


NOTE:  More  than  one  card  may  be  read  here. 


FIELD 


REMARKS 


CONTENTS 

Global  location  of  Ith  approximating  variable. 


•  If  INXLOC  *  0  in  BLOCK  J,  this  data  is  not  read.  In  this  case,  the 
data  is  defaulted  to  be  the  global  locations  of  the  design  variables 
( IDSGN  values  in  BLOCK  G). 


$ 


5wv.  .lmwISvsKw 


DATA  BLOCK  H  OMIT  IF  NXAPRX  =  0  IN  BLOCK  B  OR  INFLOC  -  0  IN  BLOCK  J 

DESCRIPTION:  Global  locations  of  functions  to  be  approximated. 


FORMAT  AND  EXAMPLE 

1  2  3  A  5  6  7  8  FORMAT 

LOCFl  LOCF2  LOCF3  LOCF4  .  8110 

3  5  6  4 


NOTE:  More  than  one  card  may  be  read  here. 

FIELD  CONTENTS 

1-8  LOCI:  Global  location  of  Ith  function  to  be  approximated. 

REMARKS 

•  If  INXLOC  •  0  in  BLOCK  J,  this  data  is  not  read.  In  this  case,  the 
data  is  defaulted  to  be  the  global  locations  of  the  objective  function 
(IOBJ  in  BLOCK  E)  followed  by  the  global  locations  of  the  constrained 
parameters  (ICON,  JCON  in  BLOCK  I). 


DATA  BLOCK  _N  OMIT  IF  NXS  =  0  IN  BLOCK  J 
DESCRIPTION:  X-vectors  for  approximate  optimization. 


FORMAT  AND 

1 

EXAMPLE 

2 

3 

4  5 

6 

7 

8 

FORMAT 

XII 

XI2 

XI3 

XI4 

. . . 

•  •  • 

. . . 

8F10 

4. 

NOTE :  NXS 

15. 

Sets  of 

data 

are  read  here. 

NOTE:  More  than  one  card  may  be  read  for  each  set  of  data. 

FIELD  CONTENTS 

1-8  XI J:  Jth  value  of  Ith  X-vector,  J  =  1, NXAPRX. 


DATA  BLOCK  0_  OMIT  IF  NXFS  -  0  IN  BLOCK  J 

DESCRIPTION:  X~F  pairs  of  information  for  approximate  optimization. 
FORMAT  AND  EXAMPLE 


1 

2  3 

4 

5 

6  7  8 

FORMAT 

XI 

X2  X3 

X4 

.  . 

...  ...  ... 

8F10 

2. 

18. 

Y1 

Y2  Y3 

Y4 

Y5 

•  •  •  »*.  ... 

7200. 

416.667  .914495 

18418.419 

NOTE: 

NSFS 

sets  of  data  are 

read  here. 

NOTE: 

More 

than  one  card  may 

be  required 

for 

XI  or  YI. 

NOTE: 

NXAPRX  values  of  X  and 

NF  values  of 

Y 

are  read  for  each 

set  of  data 

NOTE: 

Many 

significant  digits  are  desirable. 

FIELD 

CONTENTS 

1-8 

XI : 

Ith  value  of 

X,  I  -  1, NXAPRX. 

1-8 

YI: 

Ith  value  of 

Y,  I  -  1 ,NF. 
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DATA  BLOCK 


P  OMIT  IF  NSV  ’  fl  IN  BLOCK  B 


DESCRIPTION:  Sensitivity  objectives  (function  values). 
FORMAT  AND  EXAMPLE 


FORMAT 


NSOBJ  IPSENS 


NSN3  NSN4  NSN5 


Two  or  more  cards  are  read  here. 


FIELD 


CONTENTS 


NSOBJ:  Number  of  separate  objective  functions  to  be  calculated  as 

functions  of  the  sensitivity  variables. 

IPSENS:  Print  control.  If  IPSENS. GT.O,  detailed  print  will  be  called 
at  each  step  in  the  sensitivity  analysis.  DEFAULT  *  No 
print. 

NSNI:  Global  variable  number  associated  with  the  sensitivity 

objective  functions. 


REMARKS 


•  More  than  eight  sensitivity  objectives  are  allowed.  Add  data  cards  as 
required  to  contain  data. 


I 


I 

$ 


DATA  BLOCK  £  OMIT  IF  NSV  -  0  IN  BLOCK  B 
DESCRIPTION:  Sensitivity  variables. 

FORMAT  AND  EXAMPLE 

12  3  4  5 

I SENS  NSENS 


8  FORMAT 


SNS 1  SNS2  SNS3  SNS4  .  8F10 

200.  100.  150.  250. 

NOTE:  Read  one  set  of  data  for  each  of  the  NSV  sensitivity  variables. 

NOTE:  Two  or  more  cards  are  read  for  each  set  of  data. 

FIELD  CONTENTS 

1  ISENS:  Global  variable  number  associated  with  the  sensitivity 

variable . 

2  NSENS:  Number  of  values  of  this  sensitivity  variable  to  be  read  on 

the  next  card. 

1-8  SENSI:  Values  of  the  sensitivity  variable.  I  »  l, NSENS.  I  **  1 

correspond  to  the  nominal  value. 


REMARKS 


■  More  than  eight  values  of  the  sensitivity  variable  are  allowed.  Add 
data  cards  as  required  to  contain  the  data. 


DATA  BLOCK  R_  OMIT  IF  N2VAR  =  0  IN  BLOCK  B 

DESCRIPTION:  Two  variable  function  space  control  parameters. 

FORMAT  AND  EXAMPLE 

1  2  3  A  5  FORMAT 

N2VX  M2VX  N2VY  M2VY  IP2VAR  5110 

1  A  2  5  0 

FIELD  CONTENTS 

1  N2VX:  Global  location  of  the  X-variable  in  the  two  variable 

function  space. 

2  M2VX:  Number  of  values  of  X  to  be  considered. 

3  N2VY:  Global  location  of  the  Y-variable  in  the  two  variable 

function  space. 

A  M2VY:  Number  of  values  of  Y  to  be  considered. 

5  IP2VAR:  Print  control.  If  IP2VAR.GT.0,  detailed  print  will  be  called 

at  each  step  (each  X-Y  combination).  DEFAULT  «  No  print. 


DATA  BLOCK  S_  OMIT  IF  N2VAR  =  0  IN  BLOCK  B 

DESCRIPTION:  Objective  functions  of  the  two  variable  function  space  study. 
FORMAT  AND  EXAMPLE 

1  2  3  A  5  6  7  8  FORMAT 

NZ1  NZ2  NZ3  NZA  NZ5  .  8110 

3  A  5  6  7 

FIELD  CONTENTS 

1-8  NZ1 :  Global  location  corresponding  to  the  Ith  function  of  X  and  Y 

to  be  calculated.  N2VAR  values  are  read  here. 

REMARKS 

•  More  than  eight  objective  functions  are  allowed.  Add  data  cards  as 
required  to  contain  the  data. 


DATA  BLOCK  T_  OMIT  IF  N2VAR  =  0  IN  BLOCK  B 

DESCRIPTION:  Values  of  the  X-variable  in  a  two  variable  function  space  study. 
FORMAT  AND  EXAMPLE 

12  3  4  5  6  7  8  FORMAT 

XI  X2  X3  X4  .  8F10 

0.5  1.0  1.5  2.0 

FIELD  CONTENTS 

1~8  XI:  Values  of  the  X-variable  in  the  two  variable  function  space. 

M2VX  values  are  read  here. 

REMARKS 

•  More  than  eight  values  are  allowed.  Add  data  cards  as  required  to 
contain  the  data. 


DATA  BLOCK  £  OMIT  IF  N2VAR  -  0  IN  BLOCK  B 

DESCRIPTION:  Values  of  the  Y-variable  in  a  two  variable  function  space  study. 
FORMAT  AND  EXAMPLE 

1  2  3  4  5  6  7  8  FORMAT 

Y 1  Y2  Y3  Y4  Y5  .  8F10 

4.0  8.0  12.0  16.0  20.0 

FIELD  CONTENTS 

1-8  YI:  Values  of  the  Y-variable  in  the  two  variable  function  space. 

N2VY  values  are  read  here. 

REMARKS 

•  More  than  eight  values  are  allowed.  Add  data  cards  as  required  to 
contain  the  data. 


DATA  BLOCK  V_ 

DESCRIPTION:  COPES  data  'END'  card. 
FORMAT  AND  EXAMPLE 


FORMAT 


FIELD 


CONTENTS 

The  word  'END'  in  columns  1-3 


REMARKS 


•  This  card  MUST  appear  at  the  end  of  the  COPES  data 

•  This  ends  the  COPES  input  data 

•  Data  for  the  user-supplied  routine,  ANALIZ,  follows  this. 


5.2  TRANSONIC  ANALYSIS  INPUT  DATA 


NOTE:  Excluding  literal 

cards , 

all  input  data  cards  are  7F10. 

format. 

CARD 

NUMBER 

CARD 

COLUMN 

VARIABLE 

NAME 

DESCRIPTION 

Card  1 -A 

1-80 

TTLE 

Configuration  or  run  title  to 

identify 

Card  2-A 

1-10 

CASE 

graphic  and  printed  output. 

-  1.  Isolated  Body,  used  for 

input  check  of 

complex  body  definition.  No  flow 
solution.  (Omit  cards  3~A,  4-A  and 
all  cards  -C,  -W). 

=  2.  Isolated  Wing  (omit  all  cards  ~C,  -B). 
»  3.  Wing-Body  (omit  cards  ~C). 

-  4.  Isolated  wing-canard  (omit  all  cards 

-B). 

*  5.  Wing-Body-Canard  (omit  cards  -V). 

•  6.  Wing-T-Tail. 

-  7.  Wing-Body-T-Tail. 

11-20  AMACH  Mach  Number  (AMACH  <  1.0). 

21-30  AOA  Angle-of-Attack  (degrees). 

31-40  WPO  <.  -3.  Same  as  WPO  -  -2.  plus  omit  grid  and 

body  information  output. 

<  -2.  Same  as  WPO  “  -1.  plus  omit  body 
output. 

£  -1.  Same  as  WPO  ■  0.  plus  omit  Mach 
chart  output. 

£  0.  No  crude  grid  output. 

_>  1.  Crude  grid  output  for  diagnostic 
purposes . 

2.  Same  as  WPO  ”  1.  plus  boundary  layer 
information. 

_>  3.  Same  as  WPO  ■  2.  plus  print  C  values 
off  wing  and  canard. 

_>  4.  Same  as  WPO  “  3.  plus  print  fine  grid 
boundary  conditions. 


CARD 

NUMBER 


CARD 

COLUMN 


VARIABLE 

NAME 


DESCRIPTION 


Card  2-A  >  5.  Same  as  WPO  ■  4.  plus  print  solution 

(contd)  convergence  information  for  every  spanwise 

plane. 

41-50  AXIT  Number  of  initial  crude  grid  iterations. 

51-60  AXITF  Number  of  crude/fine  grid  iteration  cycles. 

61-70  VISMOD  *  1.  No  viscous  effects. 

•  2.  Viscous  effects  computed  at  end  of 
inviscid  analysis. 

=  3.  Inviscid/viscous  interaction. 

Card  3-A  1-10  FSAVE  *  1.  Save  the  flow  solution  on  Unit  99. 

t  1.  Do  not  save  the  flow  solution. 

11-20  FSTRT  -  1.  Restart  the  flow  solution  from 

Unit  98. 

t  1.  Do  not  restart  the  flow  solution. 

NOTE:  See  Usage  Note  1. 

21-30  CNVTST  Convergence  test  based  on  average  flow 

solution  correction  (CCAV).  Default  is 
1.0E-06. 

31-40  FCASM2  -  1.  Construct  wing-canard  or  wing-tail 

grid,  but  run  canard-off  or  tail-off 
solution.  Mill  read  canard  or  tail 
input  (CASE  -  4.  through  7.  on  Card 
2-A),  construct  grids  and  then  do  wing 
alone  or  wing-body  solution. 

£  1.  Do  not  "turn  off"  canard  or  tail. 

41-50  PCTLE  Crude  grid  leading  edge  spacing  tolerance. 

Default  is  0.01,  minimum  allowed  is  0.005. 
NOTE:  See  Usage  Note  2. 

51-60  FCGRD  “  1.  Places  the  canard  surface  in  the 

vertical  mid-location  of  the  canard 
fine  grid  system.  Appropriate  for 
lightly  or  negatively  loaded  canards. 


CARD 

NUMBER 

Card  3-A 
(contd) 


Card  4-A 


NOTE:  Omit 

Card  1-C 


CARD 

COLUMN 

VARIABLE 

NAME 

DESCRIPTION 

61-70 

FTWAKE 

t  1.  Places  the  canard  surface  at  the  lower 

quarter  vertical  location  in  the 

canard  fine  grid  system. 

=  0.  No  wake  deflection  (FTWAKE  defaults  to 

1-10 

SREF 

10000). 

1.  Wake  deflection  updated  every  FTWAKE 

iterations.  Recommended  value  is  8. 

Reference  area,  if  0.,  code  will  calculate. 

11-20 

AMAC 

Mean  aerodynamic  chord,  if  0.,  code  will 

21-30 

ALAM 

calculate . 

Reference  taper  ratio,  if  0.,  code  will 

31-40 

XMOM 

calculate. 

X-position  for  pitching  moment  reference. 

41-50 

ZMOM 

Z-position  for  pitching  moment  reference. 

51-60 

RE 

Reynolds  Number  x  10  based  on  AMAC. 

61-70 

FYINT 

DEFAULT  -  10. 

«  0.  Nondimens ional  ordinate  spanwise 

Card  1-C 

for  CASE  <  4 

interpolation. 

t  1.  Physical  ordinate  spanwise 

interpolation. 

NOTE:  See  Usage  Note  3. 

• 

1-10 

ASECT 

Number  of  streamwise  sections  defining  the 

11-20 

ANIN 

canard  planform  (2.  ^  ASECT  20.).  ASECT 

for  canard  plus  ASECT  for  wing  must  not 

exceed  20. 

Number  of  ordinates  defining  each  canard 

21-30 

ANOSW 

section  (ANIN  <  60). 

»  0.  Sharp  nose  canard  section. 

31-40 

ZWINGC 

-  1.  Blunt  nose  canard  section. 

Z-position  of  canard  (waterline). 

31-40 


CARD 

NUMBER 


CARD 

COLUMN 


VARIABLE 

NAME 


DESCRIPTION 


Card 

1-C 

41-50 

XTRNC 

Transition  location  for  canard,  streamwise. 

(con 

td) 

>  0.  fraction  of  chord. 

<  0.  physical  distance  from  leading  edge. 

•  0.  default  to  fixed  chord  fraction  of 

0.05. 

NOTE:  See  Usage  Note  4. 

51-60 

CINSDS 

Canard  incidence,  degrees,  incorporated  as 

added  twist  at  each  input  station. 

61-70 

SHIFTC 

Grid  shift  for  CASE  *  4  or  CASE  -  5. 

NOTE:  See  Usage  Note  5. 

Card 

l-w 

1-10 

ASECT 

Number  of  streamwise  sections  defining  the 
wing  planform  (2.  <_  ASECT  <_20.). 

ASECT  for  canard  plus  ASECT  for  wing  must 

not  exceed  20. 

11-20 

ANIN 

Number  of  ordinates  defining  each  wing 

section  (ANIN  60.). 

21-30 

ANOSW 

-  0.  Sharp  nose  wing  sections. 

-  1.  Blunt  nose  wing  sections. 

31-40 

ZWING 

Z-position  of  wing  (waterline). 

41-50 

XTRNW 

Transition  location  for  wing,  streamwise. 

Same  usage  as  for  XTRNC  (See  Card  1-C  and 

Usage  Note  4.). 

51-60 

WINSDS 

Wing  incidence,  degrees,  incorporated  as 

added  twist  at  each  input  station. 

61-70 

SHIFTW 

Grid  shift  for  CASE  <  4. 

NOTE:  See  Usage  Note  5. 

NOTE : 

Omit 

card  1-V 

for  CASE  < 

6. 

Card 

I-V 

1-10 

ASECT 

Number  of  streamwise  sections  defining  the 

vertical  stabilizer  (2.  <_  ASECT  <_  5 .  ) . 

11-20 

ANIN 

Number  of  ordinates  defining  each  input 

section  (ANIN  <  60.). 
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CARD 

NUMBER 


CARD 

COLUMN 


VARIABLE 

NAME 


DESCRIPTION 


Card 

l-v 

21-30 

ANOSW 

■  0.  Sharp  nose  sections 

(contd) 

“  1.  Blunt  nose  sections. 

31-40 

SHIFTT 

Grid  shift  for  CASE  >  5. 

NOTE:  See  Usage  Note  5. 

NOTE: 

Omit 

Card  Set 

2-C 

through  5-C  for  CASE  <  4.  Card  Set  2~C  through  5~C 

is  repeated  ASECT 

(for  canard)  times. 

Card 

2-C 

1-10 

XPL 

Canard  section  leading  edge  (X-value). 

11-20 

YP 

Canard  section  span  position  (Y-value). 
First  Y-Value  must  be  0.0  (symmetry  plane), 

even  for  wing-body  case. 

21-30 

XPT 

Canard  section  trailing  edge  (X-value). 

31-40 

TWIST 

Canard  section  local  incidence  (twist 

angle),  degrees,  added  to  input  section. 

41-50 

AKODE 

*  0.  Section  ordinates  identical  to 

preceding  section  (omit  cards  3~C 
through  5-C). 

-  1.  New  section  definition  expected  on 

cards  4-C  and  5-C. 

Card 

3-C 

1-70 

XINWC 

Canard  section  x/c  coordinates  (cards  3-C 

defined  only  for  first  canard  section,  ANIN 

values) . 

Card 

4-C 

1-70 

YINU 

Canard  section  upper  surface  y/c 

coordinates  (ANIN  values). 

Card 

5-C 

1-70 

YINL 

Canard  section  lower  surface  y/c 

coordinates  (ANIN  values). 

NOTE: 

Card 

Set  2-W 

through  5-W 

is  repeated  ASECT  (for  wing)  times. 

Card 

2-W 

1-10 

XPL 

Wing  section  leading  edge  (X  value). 

11-20 

YP 

Wing  section  span  position  (Y-value). 

First  Y-value  must  be  0.0  (symmetry  plane), 

even  for  wing-body  case. 

21-30 

XPT 

Wing  section  trailing  edge  (X-value). 

31-40 

TWIST 

Wing  section  local  incidence  (twist  angle), 
degrees,  added  to  input  section. 

CARD 

NUMBER 


CARD 

COLUMN 


VARIABLE 

NAME 


DESCRIPTION 


Card  2-W  41-50  AKODE  “  0.  Section  ordinates  identical  to 

(contd)  preceding  section  (omit  cards  3~W 

through  5-W). 

-  1.  New  section  definition  expected  on 
cards  4-W  and  5-W. 


Card 

3-W 

1-70 

XINW 

Wing  section  x/c  coordinates  (cards  3-W 

defined  only  for  first  wing  section,  ANIN 

values  expected). 

Card 

4-W 

1-70 

YINU 

Wing  section  upper  surface  y/c  coordinates 

(ANIN  values). 

Card 

5-W 

1-70 

YINL 

Wing  section  lower  surface  y/c  coordinates 

(ANIN  values). 

NOTE:  Omit  Card  Set  2-V  through  4-V  for  Case  <  6.  Card  Set  2-V  through  4-V 

is  repeated  ASECT  (for  vertical  stabilizer)  times. 


Card  2~V 

1-10 

XPLV 

Vertical  stabilizer  section  leading  edge 

(X-value) . 

11-20 

ZPV 

Vertical  stabilizer  section  Waterline 

(2-value). 

21-30 

XPTV 

Vertical  stabilizer  section  trailing  edge 

(X-value) . 

31-40 

AKODE 

“  0.  Section  ordinates  identical  to 

preceding  section  (omit  cards  3-V  and 

4-V). 

=  1.  New  section  definition  expected  on 

card  4-V. 

Card(s)  3-V 

1-70 

XINV 

Vertical  stabilizer  section  x/c  coordinates 

(Card(s)  3-V  defined  only  for  first  input 
section,  ANIN  values). 

Card(s)  4-V 

1-70 

YINV 

Vertical  stabilizer  y/c  coordinates  (ANIN 
values).  Code  allows  symmetrical  sections 

CARD  CARD  VARIABLE 

NUMBER  COLUMN  NAME  DESCRIPTION 

NOTE:  Omit  card  set  1-B  through  13-B  for  CASE  ”  2,  CASE  ■  4,  or  CASE  "  6. 

Card  1-B  1-10  BKOD  -  1.  Infinite  cylinder  (only  RADIUS  need  be 

input) . 

-  1.  Same  as  BKOD  -  1.  No  embedded  body 

grid.  Crude  grid  body  representation 
only. 

-  2.  Simple  axisymmetric  body  definition 

requested  (input  XBIN,  RIN  on  card(s) 
2-B  and  3-B. 

—2.  Same  as  BKOD  "2.  No  embedded  body 
grid.  Crude  grid  body  representation 
only. 

-  3.  Complex  body  definition  requested 

(input  Quick-Geometry  model  on  card(s) 
4-B  through  13-B).  Detailed  model 
interrogation  for  checking  input. 

“~3.  Same  as  BKOD  “3.  No  embedded  body 
grid.  Crude  grid  body  representation 
only. 

NOTE:  BKOD  >  0.  is  available  for  body  input  checkout  only.  Flow  solution  is 

available  for  BKOD  <  0.  only. 

NOTE:  See  Usage  Note  6. 

11-20  BNOSE  Body  nose  (X-value)  used  for  BKOD  t  1. 

21-30  BTAIL  Body  tail  (X-value)  used  for  BKOD  *  1. 

Card  1-B  31-40  BNIN  Number  of  axisymmetric  body  coordinates  to 

be  input.  BNIN  <  60.  (for  BKOD  -  ±2. 


■ 

CARD 

CARD 

VARIABLE 

V 

NUMBER 

COLUMN 

NAME 

DESCRIPTION 

j'l 

1 

Card  1~B 

41-50 

RADIUS 

Cylinder  radius  for  BKOD  -  ±1.  only. 

8 

g 

(contd) 

51-60 

ANOSB 

~  0.  Sharp  nose  body. 

=  1.  Blunt  nose  body. 

S 

Used  for  BKOD  =  ±2.  only. 

I 

NOTE:  Omit 

card  sets 

2-B  and  3-B 

for  BKOD  -  ±1.  or  BKOD  -  ±3. 

Card(s)  2-B 

1-70 

XINB 

Axisymmetr ic  body  X-coordinates  (BNIN 

values ) . 

i 

Card(s)  3-B 

1-70 

RIN 

Axisymmetr ic  body  radii  (BNIN  values). 

| 

NOTE:  Omit 

card  sets 

4-B  through 

13-B  for  BKOD  »  ±1.  or  BKOD  -  ±2. 

Card  4-B 

1-70 

VTITLE 

Quick-Geometry  model  title. 

5 

Card  5-B 

1-10 

ACSM 

Number  of  distinct  cross-section  models 

a 

(ACSM  card  sets  6-B  and  7-B  will  follow). 

J 

Card  6-B 

1-10 

ADUM 

Running  count  of  current  cross-section 

model  (1-ACSM). 

a 

11-20 

AARC 

Number  of  arcs  in  current  cross-section 

model  (AARC  Card(s)  7-B  will  follow). 

K 

| 

21-60 

CTITLE 

Title  or  descriptor  of  current  cross- 

section  model. 

P 

Card  7-B 

1-8 

ARCNAM 

Arc  or  component  name. 

h 

P 

11-14 

ASHAPE 

Arc  or  component  shape. 

» 

21-28 

PNTNAM(l) 

Control  point  name  for  beginning  of  this 

I 

arc. 

i 

31-38 

PNTNAM(2) 

Control  point  name  for  termination  of  this 

| 

arc. 

j 

41-48 

PNTNAM(3) 

Slope  control  point  name  for  this  arc,  if 

required. 

» 

m 

Card  8-B 

1-10 

ANTCSM 

Number  of  cross-section  models  to  define 

entire  body  (ANTCSM  card(s)  9-B  will 

follow) . 

\ 
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CARD 

NUMBER 

CARD 

COLUMN 

VARIABLE 

NAME 

Card  9-B 

1-10 

ADUM 

11-20 

AMODEL 

21-30 

XCSMS1 

31-40 

XCSMS2 

Card  10-B 

1-10 

BLINE 

11-20 

ALIAS 

Card 

11-B 

1-10 

BLSEG 

11 

BYORZ 

12-19 

BNAME 

Card 

12-B 

1-4 

S SHAPE 

11-20 

D(  1 ) 

21-30 

D(2) 

31-40 

D(3) 

►— » 

1 

LA 

O 

D(4) 

51-60 

D(  5  ) 

61-70 

D(6) 

Card 

13-B 

11 

BYORZ 

12-19 

BNAME 

21 

AYORZ 

22-29 

ANAME 

DESCRIPTION 

Running  count  of  current  cross-section 
model  (1-ANTCSM). 

Index  corresponding  to  already  defined 
cross-section  models  (between  1  and  ACSM). 
Starting  X-station  for  current 
cross-section  model. 

Ending  X-station  for  current  cross-section 
model . 

Number  of  body  line  models  to  be  defined  by 
segments  (BLINE  card  set  11-B  and  12-B 
follow) . 

Number  of  body  line  models  to  be  aliased. 
(Input  ALIAS  card(s)  13-B  below). 

Number  of  segment(s)  defining  body  line 
mode  1 . 

The  letter  Y  or  Z  indicates  which  data 
definition  is  to  follow. 

Body  line  name  to  be  defined. 

Segment  shape. 

X-station  for  beginning  of  segment. 

Y  or  Z  value  corresponding  to  D(l). 
X-station  for  termination  of  segment. 

Y  or  Z  value  corresponding  to  D(3). 
X-station  for  segment  slope  control  point. 

Y  or  Z  value  corresponding  to  D(5). 

The  letter  Y  or  Z  indicates  which  data 
definition  is  to  follow. 

Body  line  name  to  be  defined. 

The  letter  Y  or  Z  indicates  which 
definition  is  to  be  used  for  aliasing. 

Body  line  name  to  which  BNAME  is  aliased. 
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5.3  USAGE  NOTES 

Note  1.  Saved  Solution 

When  an  analysis  consists  of  both  crude  only  and  then  crude/fine 
iterations,  a  saved  solution  is  written  on  unit  99  immediately  before  the 
crude/  fine  iterations  begin.  At  the  end  of  the  crude/fine  iterations, 
the  previous  saved  solution  is  overwritten  with  the  most  recent  results. 
In  this  way,  an  abnormal  termination  of  the  crude/fine  iterations  will 
have  the  crude  iteration  results  saved  on  unit  99.  With  an  abnormal 
termination,  appropriate  job  control  cards  may  be  needed  to  "permanently" 
save  the  unit  99  data. 

Note  2.  PCTLE  Parameter 

The  code  will  shift  the  input  geometry  to  find  a  streamwise  location 
in  the  crude  grid  where  the  first  points  at  each  span  station  on  the 
wing  and  canard  is  not  "too  close"  to  the  leading  edge.  (Riegel's  rule 
is  not  used  for  the  boundary  conditions.)  The  PCTLE  parameter  is  the  re¬ 
quired  minimum  distance  from  the  leading  edge.  See  also  Usage  Note  5, 
Grid  Shift  Parameters. 

Note  3.  FYINT  Parameter 

The  input  wing  and  canard  section  ordinates  (YINU  AND  YINL)  are  lin¬ 
early  interpolated  from  the  input  span  stations  to  the  analysis  span  sta¬ 
tions.  The  method  of  interpolation  is  controlled  by  the  parameter  FYINT. 
If  FYINT  *  0.  nondimens ional  spanwise  interpolation  of  the  ordinates  is 
used : 

(z/c)y  -  (z/c)yl(l.-R)  +  <z/c)y2(R) 

If  FYINT  *  1.  physical  spanwise  interpolation  of  the  ordinates  is  used: 

(z/c)y  =  (z/c)yl(l.-R)(cyl/cy)  +  (z/c)y2(R)(cy2/cy) 
where 

c  ■  local  chord  Y  ■  interpolated  span  station 

R  =  (Y-Y1)/(Y2-Y1)  Y1,Y2  «  input  span  stations 

z/c  =  nondimens ional  section  ordinate 


The  first  formula  is  the  more  usual  analysis  code  interpolation,  while 
the  second  formula  corresponds  to  manufacturing  lofting  methods.  The 
difference  between  the  two  formulas  becomes  greater  for  more  highly 
tapered  wings  and  fewer  input  span  stations. 

Note  4.  Boundary  Layer  Transition 

The  fraction  of  chord  designation  results  in  transition  specified  at 
the  same  percentage  chord  location  at  ail  span  stations.  Thus,  the  phys¬ 
ical  distance  from  the  leading  edge  decreases  as  the  local  chord  length 
decreases.  The  physical  distance  designation  results  in  transition  spec¬ 
ified  at  the  same  fixed  distance  from  the  leading  edge.  Thus,  the 
percentage  chord  from  the  leading  edge  increases  as  the  local  chord 
length  decreases. 

Note  5.  Grid  Shift  Parameters 

Variables  SHIFTC,  SHIFTW  and  SHIFTT  can  be  used  to  set  the  initial 
streamwise  placement  of  the  configuration.  (The  grid  is  fixed  in  loca¬ 
tion  and  the  input  geometry  is  shifted.)  The  code  will  determinate  the 
appropriate  shift  with  a  search  procedure.  After  the  first  analysis  run, 
the  resulting  shift  parameter  can  be  input  to  eliminate  the  search  in 
subsequent  runs.  The  shift  parameter  can  also  be  used  to  evaluate  dif¬ 
ferent  grid  placements.  If  the  grid  placement  search  should  fail,  the 
shift  parameter  can  be  used  to  start  the  search  at  different  locations, 
possibly  finding  a  satisfactory  grid  placement. 

Note  6.  Body/Fuselage  Geometry  Model 

The  present  method  allows  complex  three-dimensional  geometries  to  be 
input,  processed  and  converted  into  a  suitable  array  of  boundary  condi¬ 
tions  for  analysis.  Although  the  input  or  modeling  of  complex  body 
shapes  is  extremely  error  prone  and  certain  applications  might  not  war¬ 
rant  this  level  of  effort,  it  is  necessary  in  aircraft  application  when 
fuselage  contours  (e.g.,  canopies,  fairings)  are  required.  The  code  does 
not  provide  an  embedded  fine  grid  body  analysis.  The  embedded  body  grid 
option  is  provided  to  allow  a  detailed  body  input  definition  checkout. 


This  will  be  most  useful  if  the  graphical  output  is  available.  The  fol¬ 
lowing  discussion  is  excerpted  from  Ref  2.  Additional  information  and 
examples  will  be  found  in  the  reference. 

The  complex  fuselage  modeling  system  has  been  named  "Quick-Geometry" 
by  its  developers,  Vachris  and  Yaeger  (Ref  3).  A  detailed  User's  Guide 
for  the  Quick-Geometry  system  can  be  found  in  the  appendix  of  Ref  14. 

This  system  was  originally  developed  for  the  geometric  modeling  of  wing- 
body  shapes.  Since  only  fuselage  shapes  are  of  concern  here,  many  of  the 
more  sophisticated  options  including  fillets  and  patches  will  not  be  de¬ 
scribed  in  the  paragraphs  which  follow.  In  addition,  if  Ref  3  and  14  are 
being  used  to  augment  the  modeling  description  provided  herein,  it  should 
be  noted  that  the  input  format  has  been  modified  to  be  more  consistent 
with  that  of  the  basic  transonic  wing-body  code. 

The  geometry  package  requires  that  certain  body  line  and  cross-section 
lines  be  defined.  The  body  lines  and  cross-section  lines  may  be  likened 
to  the  stringers  and  bulkheads,  respectively,  used  in  fuselage  construc¬ 
tion.  These  line  models  are  defined  by  a  combination  of  simple  curves 
(i.e.,  lines,  ellipses,  cubics).  They  are  taken  together  to  provide  a 
continuous  analytical  model  of  the  surface  geometry.  Slopes  and  normals 
are  developed  analytically.  Either  discontinuous  intersections  or  smooth 
fairings  can  be  modeled  and  enforced. 

Two  different  coordinate  systems  are  employed.  Geometry  definition  is 
performed  in  a  Cartesian  coordinate  system  (X,Y,Z),  while  interrogation 
of  the  model  for  body  boundary  conditions  is  performed  in  cylindrical 
coordinates  (X,R,0).  This  results  in  the  use  of  a  plane  of  symmetry  map 
axis,  the  height  of  which  usually  corresponds  to  the  position  of  the 
max-half-breadth  line.  It  is  required  that  the  configuration  radius  at 
any  cross-sectional  cut  be  a  single  valued  function  of  the  angle  0. 

These  definition  lines  and  coordinate  systems  are  illustrated  in  Figure 
14. 

Marconi,  F.  and  Yaeger,  L. ,  "Development  of  A  Computer  Code  for  Calcu¬ 
lating  the  Steady  Super/Hypersonic  Inviscid  Flow  Around  Real  Configura¬ 
tions,"  NASA  CR-2676  (Vol.  II),  May  1976. 
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FiBure  14.  Quick  geometry  model  lines  and  coordinate  systems  (from  Ref  2). 


A  minimum  of  four  body  lines  are  required  for  the  simplest  fuselage. 
These  are:  top  centerline,  bottom  centerline,  max-half-breadth  line,  and 
the  map  axis.  Each  body  line  must  be  defined  by  both  its  Y  and  Z  values 
over  the  full  range  of  X  (between  fuselage  nose  and  tail).  Similarly,  a 
minimum  of  two  cross-section  line  segments  are  required  for  each  differ¬ 
ent  cross-section  line  model.  These  are  body  upper,  and  body  lower. 

Both  body  lines  and  cross-section  lines  are  specified  by  defining  key 
arc  or  segment  shapes  and  their  accompanying  limiters.  The  segment  shape 
boundary  conditions  used  to  determine  the  coefficients  of  the  slope  equa¬ 
tion  are  the  origin  point,  termination  point,  and  slope  control  point. 

The  slope  control  point  lies  at  the  intersection  of  the  line  which  is 
tangent  to  the  segment  shape  at  the  origin  point  and  the  line  which  is 
tangent  to  the  segment  shape  at  the  termination  point  (see  Figure  15). 

The  slope  control  point  is  a  very  convenient  way  of  specifying  slope  con¬ 
ditions.  In  particular,  it  allows  for  the  simultaneous  specif ication  of 
slope  conditions  at  both  ends  of  the  segment. 


am asKPsassa 


Figure  15  is  a  schematic  illustrating  the  component  build-up  of  a  par¬ 
ticular  body  line  and  cross-section  line  model.  Naturally,  LINE  segments 
do  not  require  a  slope  control  point.  In  this  case,  the  portion  of  the 
body  top  centerline  illustrated  requires  four  body  line  segments  and  the 
cross-section  is  constructed  with  two  arcs  (two  is  the  minimum  number 
a  1  lowed) . 

Components  of  body 
(top  center)  line  model 


Canopy  slope  control  point 


Figure  15.  Quick  geometry  body  and  crots -section  line  models  (from  Ref  2). 

The  arc  shapes  used  for  defining  a  cross-section  line  model  are  listed 
in  Table  1.  They  are  input  in  an  order  which  starts  at  the  body  bottom 
centerline  and  proceeds  to  the  body  top  center  line.  The  segment  shapes 
used  for  defining  a  body  line  model  are  listed  in  Table  2. 

Cross-section  arcs  are  input  in  their  order  of  appearance.  However, 
body  line  segments  are  defined  along  with  an  index  which  establishes 
their  order  in  the  X  direction.  In  addition,  body  lines  may  be  aliased 
to  other  body  lines.  This  allows  two  body  line  models  to  have  identical 
mathematical  representations  without  repeating  the  body  line  segment 
input.  For  example,  the  Z  value  of  the  map  axis  (ZMAP)  is  typically  the 
Z  value  of  the  maximum  half-breadth  (ZMHB).  The  two  are  made  identical 
by  aliasing  the  two  body  line  model  names  ZMAP  and  ZMHB. 


TABLE  1.  CROSS-SECTION  ARC  SHAPES 


SHAPE 

KEYWORD 

EQUATION 

LINE 

LINE 

Ay  +  Bz  +  C  -  0 

ellipse 

(CONCAVE  TO  ORIGIN) 

ELLI 

<y-y0>2  (z-zQ)2 

-  +  —  -  0 

A2  B2 

ELLIPSE 

(CONVEX  TO  ORIGIN) 

ELLO 

SAME  AS  ELLI. 

TABLE  2.  BODY  LINE  SEGMENTS 


SHAPE  KEYWORD 


LINE  LINE 

X  PARABOLA  XPAR 

Y PARABOLA  YPAR 

X  ELLIPSE  ELLX 

Y-ELLIPSE  ELLY 

CUBIC  CUBI 


EQUATION 


Ax  +  By  «  0 

2 

Ax  +  By  +  y  ”  0 

Ax  +  By  +  x2  -  0 

Ax  +  By  +  Cx2  +  y2  -  0 
7  7 

Ax  +  By  +  Cv  +x  «0 
Ax  +  By  +  Cx2  +  x3  -  0 


It  should  be  noted  that  cross-sections  are  defined  only  in  terms  of 
named  component  arcs  (arc  shape  table)  and  named  control  points.  On  the 
other  hand,  body  lines  are  defined  mathematically  by  coordinates  over  the 
length  of  the  configuration  for  which  they  are  required.  At  a  given  X 
station,  the  body  lines  are  interrogated  to  give  the  key  control  points 
required  to  construct  the  cross-sectional  arcs. 
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6  -  CODE  OPERATION 


The  program  reads  data  from  unit  5  and  writes  output  on  unit  6.  Units  20 
and  AO  are  used  as  scratch  files  in  the  COPES  routines.  The  scratch  file  num¬ 
bers  may  be  changed  by  changing  two  cards  at  the  beginning  of  the  COPES  pro¬ 
gram.  The  analysis  code  stores  data  on  unit  numbers  1,  A,  8,  15,  85,  90,  9A, 
98  and  99.  Unit  15  is  used  by  the  Ames  plotting  software  for  graphical  out¬ 
put.  It  has  been  left  available  to  facilitate  incorporation  of  graphic  capa¬ 
bility  by  other  users. 

The  computer  code  is  written  in  FORTRAN,  employing  the  CRAY  segmentation 
loader.  Storage  requirements  on  a  CRAY  X-MP  are  806,000  words  of  memory. 
Typical  CPU  time  for  a  wing-body-canard  analysis  with  viscous  effects  is  four 
minutes.  The  CPU  time  for  an  optimization  run  will  vary  according  to  the  com¬ 
plexity  of  the  optimization  problem,  as  discussed  in  Section  A,  COPES/COMIN 
OPTIMIZATION. 

The  segmentation  structure  is  described  in  Table  3.  The  Job  Control  Lan¬ 
guage  ( JCL)  for  operation  on  the  Ames  CRAY  is  shown  in  Figure  16.  The  job 


TABLE  3.  SEGMENTATION  STRUCTURE 


rfSF1 


JOB , JN=SAMPLE , T=660 , MFL=8SOOOO . 

-  ACCOUNT  CARD  . 

* 

COPYF, I=9IN,0=SEGDIR. 

REWIND. DN=SEGDIR. 

* 

ACCESS . DN=CANPL . PDN=CANPL10, ID=CANTATA, OWN=RFAPXA . 
UPDATE, P=CANPL,C=COMPILE. IN, ID.ED.UM. 

* 

CFT , I =COMP I LE , OFF=P , ON=A , L=0 , OPT=FULLI FCON . 

RELEASE , DN=CANPL : COMPILE . 

» 

ACCESS , DN=CANLI B , PDN=CANLB10, ID=CANTATA , OWN=RFAPXA . 
SECLDR , DW=72 , I=SEGDIR . 

RELEASE , DN=9BLD : SEGDIR . 

* 

CANTATA. 

/EOF 

PRESET=INDEF 
MAP=PART 
REDEF=IGNORE 
ABS =CANTATA 
LIB=CANLIB 
TREE 
A(B.G) 

B  (C,  D,  E) 

D(H.  I, J) 

ENDTREE 

SEGMENT=A; SAVE=ON 
MODULES=COPES 
ENDSEG 

SEGMENT=B ; SAVE=ON 
MODULES=ANALIZ 
ENDSEG 
SECMENT=C 

MODULE S=OVL 10, MLS IN 
ENDSEG 

SEGMENT=D ; SAVE=ON 
MODULE S=MLSNLS 
ENDSEG 

SEGMENT=E ; SAVE=ON 
MODULES=PLTSEG 
ENDSEG 
SEGMENT =G 

MODULES=OPT, CONMIN 
ENDSEG 
SEGMENT=H 
MODULES=STRIPK 
ENDSEG 
SECMENT=I 
MODULES=LIDRAG 
ENDSEG 

SEGMENT=J ; SAVE=ON 
MODULES=WAKEC 
ENDSEG 
/EOF 

‘COMPILE  COPES, BLKDAT 

.  ADDITIONAL  UPDATES  AS  APPROPRIATE 

/EOF 

.  INPUT  DATA  - 


Figure  16.  Job  Control  Language. 
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card  and  account  card  will  vary  with  the  individual  user.  To  save  CPU  time 
for  compilation,  a  binary  file  of  the  code  is  stored  as  a  permanent  file. 

Thus,  the  JCL  uses  an  UPDATE  command  followed  by  CFT  to  change,  compile  and 
replace  only  those  routines  being  modified  for  a  particular  run.  The  numer¬ 
ical  optimization  will  require  changes  to  the  subroutines  OGEOM,  MOD  and 
ANALIZ.  For  simple  analysis,  no  changes  need  be  made,  except  perhaps  to  BLOCK 
DATA  to  change  some  defaults.  The  UPDATE  step  always  requires  some  input. 

The  ^COMPILE  COPES, BLKDAT  card  in  Figure  16  is  shown  as  an  example  that  sat¬ 
isfies  this  requirement. 

Not  shown  in  the  JCL  are  the  cards  for  saving  and  restarting  the  flow 
solution  with  FORTRAN  files  98  and  99.  The  Applied  Computational  Fluids 
Branch  at  Ames  can  provide  the  necessary  authorization  and  information. 

If  the  code  is  converted  for  use  on  IBM  type  computers,  several  Quick- 
Geometry  variables  need  special  treatment.  These  variables  and  the  sub¬ 
routines  in  which  they  appear  are  listed  in  Table  4.  The  REAL*8  declaration 
is  needed  to  have  sufficient  word  length  for  line  model  labeling  and  increase 
the  accuracy  of  the  body  model  calculations.  The  CRAY  compiler  option  OFF-P 
is  used  to  allow  the  REAL*8  declarations  to  remain  in  the  code  and  have  no 
effect.  CRAY  compiler  directives  are  used  in  several  places.  These  begin 


TABLE  4.  REAL  *8  VARIABLES  FOR  IBM  USAGE  OF  QUICK-GEOMETRY 


SUBROUTINE  NAME 

VARIABLE  NAME 

CURVES 

A,  B,  C.  Y.  T,  X,  FACT,  RFACT,  S 

MODTV 

SUM,  ONE 

VOOTV 

C.  ONE 

QWIKDE 

CPNTNM,  COMPNM 

BlMCHK 

CPNTNM.  COMPNM,  ANAME,  BNAME,  BLMNAM,  EQUNB,  EQUNA 

BLMDEF 

CPNTNM,  COMPNM,  ANAME,  BNAME,  BLMNAM,  BLANK2 

CSMCHK 

CPNTNM,  COMPNM,  ARCNAM,  PNTNAM,  BLANK2,  ZMAPNM,  ARCNM 

CSMDEF 

CPNTNM,  COMPNM,  ARCNAM,  PNTNAM,  BLANK2,  ZMAPNM,  ARCNM,  AMAPAX 

GEMOUT 

CPNTNM,  COMPNM,  ARCNM,  PNTNAM,  BLANK 2,  AMAPNM,  ARCNM 

OSETUP 

ALABLE,  BLABLE 
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with  a  letter  C  in  the  first  column,  so  that  they  will  be  interpreted  as  com¬ 
ment  cards  on  other  computers.  Two  CRAY  library  functions,  ISMAX  and  ISAMAX 
are  used  in  the  code.  The  ISMAX  function  returns  the  index  of  the  vector  ele¬ 
ment  with  the  largest  value.  The  ISAMAX  function  returns  the  index  of  the 
vector  element  with  the  largest  absolute  value.  Their  FORTRAN  equivalents  are 
easily  incorporated.  The  overall  scalar-to-vector  CPU  time  ratio  is  2.6. 

The  program  output  listing  of  the  input  data  for  a  sample  analysis  run  is 
shown  in  Figure  17.  In  front  of  this  input  were  the  COPES  control  cards: 

TITLE 

1, 

END 

with  no  delimiters  in  between.  The  input  example  uses  NACA  0010  wing  sec¬ 
tions.  Following  the  input  data  listing  is  a  data  read  echo  showing  the  in¬ 
formation  the  program  has  read  from  the  input  cards  (Figure  18).  Portions  of 
the  grid  generation  output  are  shown  in  Figures  19  and  20. 

In  this  example,  the  flowfieid  solution  begins  with  nine  crude  grid  it¬ 
erations,  and  then  continues  with  nine  fine  grid/crude  grid  iterations.  The 
program  output  for  this  is  shown  in  Figure  21.  When  the  flow  solution  is  com¬ 
pleted,  detailed  flow  information  is  output  at  each  analysis  station  on  the 
wing  and  canard.  An  example  for  station  11  is  shown  in  Figure  22. 

Force  and  moment  output  for  the  wing,  canard  and  body  is  shown  in  Figure 
23.  Figure  24  shows  the  final  wing-body-canard  spanload  output.  All  of  the 
results  shown  in  Figure  18  through  24  were  for  the  analysis  input  data  of  Fig- 
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CANTATA  ---  CANAfId  TAIL  TRANSONIC  AERODYNAMICS 


DEVELOPED  AT  GRUMMAN  AEROSPACE  CORPORATION . BE THPAGE .  NEW  YORK 
FOR  AIR  FORCE  FLIGHT  DYNAMICS  LABORATORY.  WPAFB,  OHIO 
J.  R.  SIRBAUGH,  CONTRACT  MONITOR.  FIMM 


INPUT  DATA  LISTING 


FORWARD 

SWEPT  SAMPLE 

CASE 

5.0 

0.7 

8.0 

0.000 

9.0 

9.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

8.0 

0.0 

0.0 

0.0 

44.00 

-2.00 

0.0 

1.0 

3.0 

18  . 

1 .0 

4.00 

.05 

-2.0 

5.0 

18. 

1  .O 

-4.00 

.05 

1.0 

27  .0 

0.0 

43.0 

0.0 

1.0 

0.00 

1.25 

2.50 

5.00 

7.50 

10.0 

15.0 

20.0 

25.0 

30.0 

40.0 

50.0 

60.0 

70.0 

80.0 

90.0 

95.0 

lOO. 

0.000 

1.578 

2. 178 

2. 962 

3.500 

3.902 

4.455 

4.782 

4.952 

5.002 

4.837 

4.412 

3.803 

3.053 

2.187 

1 . 207 

0.672 

0.105 

-0.000 

-1 . 578 

-2. 178 

-2.962 

-3 . 500 

-3.902 

-4.455 

-4.782 

-4.952 

-5.002 

-4.837 

-4.412 

-3.803 

-3.053 

-2 . 187 

-1.207 

-0.672 

-0.105 

27. 

5.5 

43.0 

0.0 

0. 

33.0 

18.0 

38.0 

0.0 

0. 

54.0 

0. 

103.0 

0.0 

1.0 

0.00 

1.25 

2.50 

5.00 

7.50 

10.0 

15.0 

20.0 

25.0 

30.0 

40.0 

50.0 

60.0 

70.0 

80.0 

90.0 

95.0 

100. 

0.000 

1.578 

2.178 

2.962 

3.500 

3.902 

4.455 

4.782 

4.952 

5.002 

4.837 

4.412 

3.803 

3.053 

2 . 187 

1.207 

0.672 

0.105 

-0.000 

-1.578 

-2.178 

-2.962 

-3.500 

-3.902 

-4.455 

-4.782 

-4.952 

-5.002 

-4.837 

-4.412 

-3.803 

-3.053 

-2.187 

-1.207 

-0.672 

-0 . 105 

54.0 

5.5 

103.0 

0.0 

0.0 

57.0 

12.0 

80.00 

0.0 

0. 

58.0 

14.0 

78.0 

0.0 

0. 

38.0 

36.0 

48.0 

0.0 

0. 

-3.0 

0. 

126.0 

QUICK - 

BODY  FUSELAGE 

2  .0 

1  .0 

3.0 

AFTBOX 

ABOT 

LINE 

BDYBCL 

BDYMHB 

ASIDE 

LINE 

BDYMHB 

SIDHI 

ATOP 

LINE 

SIDHI 

BDYTCL 

2.0 

5.0 

MIDFUS 

MBOT 

LINE 

BDYBCL 

BOTCRN 

MLOCR 

ELLI 

BOTCRN 

BDYMHB 

LOSLCP 

MS  IDE 

LINE 

BDYMHB 

SIDHI 

MHICR 

ELLI 

SIDHI 

TOPCRN 

HISLCP 

MTOP 

LINE 

TOPCRN 

BDYTCL 

2.0 

1  .0 

2.0 

0.0 

27.73 

2  .0 

1.0 

27.73 

126.0 

13.0 

o.o 

Figure  17.  Analysis  code  input  data  listing  (sheet  1  of  2). 
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1.0  YCENTER 


LINE 

0.0 

0.0 

126.0 

0.0 

2.0 

ZBDYBCL 

ELLX 

0.0 

-2.2 

27.73 

-6.0 

16.75 

-6.0 

LINE 

27.73 

-6.0 

126.0 

-6.0 

2.0 

ZBDYTCL 

ELLX 

0.0 

-2.2 

27.73 

6.75 

15.12 

6.75 

LINE 

27.73 

6.75 

126.0 

6.75 

4.0 

ZBDYMHB 

LINE 

0.0 

-2 . 2 

8 . 31 

-  .6 

CUBI 

8.31 

-  .6 

16.31 

-2.0 

11.9 

0.14 

LINE 

16.31 

-2.0 

27.73 

-6.0 

LINE 

27.73 

-6.0 

126.0 

-6.0 

2. 

YBDYMHB 

CUBI 

0.0 

0.0 

27.73 

5.5 

12.7 

5.5 

LINE 

27.73 

5.5 

126.0 

5.5 

1.0 

YBOTCRN 

LINE 

0.0 

0.0 

27.73 

5.5 

1.0 

YTOPCRN 

LINE 

0.0 

0.0 

27.73 

5.5 

2.0 

YSIDHI 

CUBI 

0.0 

0.0 

27.73 

5.50 

12.7 

5.50 

LINE 

27.73 

5.5 

126.0 

5.5 

3.0 

ZSIDHI 

LINE 

0.0 

-2.2 

8.31 

0.0 

CUBI 

8.31 

0.0 

27.73 

6.75 

15.12 

4.86 

LINE 

27.73 

6.75 

126.0 

6.75 

1.0 

ZBOTCRN 

ELLX 

0.0 

-2.2 

27.73 

-6.0 

16.75 

-6.0 

1.0 

YLOSLCP 

CUBI 

0.0 

0.0 

27.73 

5.5 

12.7 

5.5 

4.0 

ZMAPAXIS 

LINE 

0.0 

-2.2 

3.4 

-1.6 

LINE 

3.4 

-1.6 

8.31 

0.0 

LINE 

8.31 

0.0 

27.73 

0.0 

LINE 
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Figure  17.  Analysis  code  input  data  lilting  (sheet  2  of  2). 
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Figure  19.  Crude  grid  generation  output  (sheet  1  of  3). 
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Figure  19.  Crude  grid  generation  output  (iheet  2  of  3). 
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Figure  19.  Crude  grid  generation  output  (sheet  3  of  3). 
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Figure  20.  Excerpt  of  fine  grid  generation  output  (sheet  1  of  2). 
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Figure  21.  Flow  solution  convergence  output  (sheet  1  of  3). 
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Figure  22.  Exemple  of  fin*  grid  solution  output  (sheet  1  of  2). 


FORWARD  SURFACE 


FORCE  AND  MOMENT  BASED  ON  SEXP=  285.847 
CL  CM  CD (PRESSURE)  CD (FRICTION) 

0.21531  0.07414  0.02890 


CL 

0.03511 

PRESSURE  DRAG  = 


CM 

0.01293 

0.00471 


CD 

0.00471 


0.00000 

BASED  ON  SREF 


FRICTION  DRAG  =  0.00000 

PITCHING  MOMENT  DUE  TO  DRAG  =  0.00084 

AFT  SURFACE 

FORCE  AND  MOMENT  BASED  ON  SEXP=  1287.256 
CL  CM  CD (PRESSURE)  CD  (FRICTION) 

0.25717  -0.09062  0.03575 

CL  CM  CD 


0.00000 

BASED  ON  SREF 


0.18885 

PRESSURE  DRAG  = 


-0.06810 

0.02626 


0.02626 


FRICTION  DRAG  =  0.00000 

PITCHINC  MOMENT  DUE  TO  DRAG  =  -0.00155 

TOTAL  FORCE  AND  MOMENT  COEFFICIENTS 
TOTAL  CL  TOTAL  CM  TOTAL  CD 

0.22396  -0.05517  0.03097 

TOTAL  PRESSURE  DRAG  =  0.03097 

TOTAL  FRICTION  DRAG  =  0.00000 

a)  Lifting  surfaces. 


BODY  LENGTH  =  126.000 

BODY  WETTED  AREA  =  5.432BE*03 

BODY  PROJECTED  AREA  =  1.2866E*03 

BODY  MAX.  CROSS-SECTIONAL  AREA  =  9.5033E-01 

X-POSITION  ABOUT  WHICH  MOMENTS  ARE  COMPUTED  =  44.000 

BODY  FORCE  AND  MOMENT  COEFFICIENTS 

BODY  CL  BODY  CM  BODY  CD 

0.03152  -0.00184  0.00465 

BODY  PRESSURE  DRAG  =  0.00465 
BODY  FRICTION  DRAG  =  0.00000 

b)  Body. 

Figure  23.  Force  and  moment  output 
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7  -  CONCLUSIONS 


The  CANTATA  code  shows  good  performance  in  the  major  areas  of  its  devel¬ 
opment.  The  embedded  grid  technique  allows  arbitrary  placement,  with  arbi¬ 
trary  overlap,  of  two  lifting  surfaces.  The  AF2YZ  algorithm  provides  a  sub¬ 
stantial  improvement  in  convergence  rate  relative  to  SLOR,  and  can  analyze 
highly  swept,  highly  tapered  fighter  wings  at  lift  coefficients  of  1.0  at 
transonic  speeds.  The  vectorization  of  the  original  scalar  oriented  code  has 
resulted  in  an  overall  scalar-to-vector  ratio  of  2.6  on  the  CRAY.  The  free 
wake  modeling  employs  a  merging  of  vortex  lines  to  treat  wake  rollup.  This 
approach  works  well  for  potential  flow  solutions  and  is  easily  transferred  to 
other  applications  of  computational  fluid  dynamics.  The  comparison  of  data 
and  analysis  indicate  that  the  free  wake  effects  are  relatively  weak,  but 
their  inclusion  generally  improves  agreement  with  test  data.  Boundary  layer 
effects  are  more  significant.  In  particular,  flow  separation  modeling  is  im¬ 
portant  for  fighter  aircraft. 


APPENDIX  A 


FREE  WAKE  CALCULATION* 

Two  modifications  were  developed  for  the  CANTATA  code  to  extend  its  capa¬ 
bilities  to  the  computation  of  "free  wake"  solutions.  The  first  modification 
introduces  a  variable  grid  spacing  that  can  handle  non-planar  wakes.  The  sec¬ 
ond  modification  introduces  an  iteration  scheme  by  which  the  flow  computation 
is  carried  out  for  several  iterations  with  a  fixed  non-planar  wake  calculated 
from  the  flowfield  of  the  previous  set  of  iterations.  Taken  together  these 
two  modifications  allow  the  computation  of  free  wake  solutions. 

A. 1  DISCRETIZATION  OF  THE  WAKE  VORTEX  SHEETS 

Completion  of  the  potential  flow  iterations  yields  the  distribution  of 
circulation  for  spanwise  grid  points  J  at  the  trailing  edge  stations  of  the 
wing  and  canard  surfaces.  For  these  distributions  the  wake  is  represented  by 
discrete  vortices  of  index  L  and  strength  T.  The  details  of  the  discretiza¬ 
tion  procedure  are  presented  in  Figure  25  for  the  wing  wake.  Computational 
grid  points  (X,Y,Z)  are  represented  by  the  indices  (l,J,K)  while  the  Lth  vor¬ 
tex  line  is  located  by  the  coordinates  (X^,  YV^,  ZVL>,  requiring  only  the 
indices  (I,L).  At  the  surface  trailing  edge  YVL  -  Yj  with  J  *  L,  i.e.,  the 
vortex  lines  originate  from  the  trailing  edge  grid  points. 

A. 2  TRACKING  OF  VORTEX  LINES 

Vortex  lines,  which  are  everywhere  tangent  to  the  local  vorticity  vector 
are  considered  the  generators  of  the  wake  vortex  sheet.  An  element  of  a  line 
vortex  must  be  tangent  to  the  vector  velocity  which  would  result  in  the  ab¬ 
sence  of  that  element.  In  addition,  a  line  vortex  element  induces  zero  veloc¬ 
ity  at  its  own  center.  In  the  limiting  case  of  a  two-dimensional  continuum  of 
line  vortices,  i.e.,  a  vortex  sheet,  these  properties  require  that  at  all 
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b)  Trailing  vortex  line  pattern 


c)  Wake  cross-section 


Figure  25.  Discretization  of  wake  vortex  sheet, 


points  the  sheet  is  tangent  to  the  mean  of  the  velocity  vectors  just  above  and 

ji. 

below  the  sheet.  If  dr  is  the  displacement  vector  of  a  line  element  corre¬ 
sponding  to  an  increment  dx  in  the  downstream  direction  we  have  for  displace¬ 
ment  between  stations  I  and  I  +  1 


a  ^ 

dr  =  VVdX/ (UV) 


(A- la) 


where 


dr  -  *  JOT»m>t  -  VV^)  *  HZVltl  L- ZVIiL) 

-ik.  Jik. 

VV“Jj( VV  +VV  )  =  iUV  +  jVV  +  k  wv 

U  1  I  »  I>  A  j  I  j  L* 


dX  =  xI+1  -  xx 


(A- lb) 
(A-lc) 
(A- Id) 


and  VV  ,  VV,  are  velocities  above  and  below  the  sheet  as  determined  from  the 
potential  solution.  Thus,  a  given  set  of  iterations  yields  values  of  the  ve¬ 
locity  vector  at  grid  points  just  above  and  below  the  assumed  wake  surface. 
These  are  linearly  interpolated  to  YV  corresponding  to  the  vortex  lines  L 

if*- 

at  station  I  to  obtain  VV.  The  equation  is  then  solved  for  (YVj+^  L  -  YV^  ^) 


and  (ZV, 


c.u  \~'T  ZV  )  and  an  updated  vortex  line  is  constructed  by  integration 

of  these  increments. 

For  the  case  of  unswept  trailing  edges  all  vortex  lines  from  a  given  sur¬ 
face  begin  at  the  same  index  I  so  that  tracking  begins  simultaneously  (the 

index  J  is  varied  for  a  fixed  station  I).  To  consider  a  case  with  a  swept 
trailing  edge  which  lies  between  the  extremes  of  I  =  ITEW_„  and  I  =  ITEUAV.  the 

MIN  MAX 

velocities  VV,  WV  are  set  to  values  which  reconstruct  the  meanline  of  the  sur¬ 
face  points  lying  between  the  two  indices.  The  tracking  process  then  begins 

at  I  =  ITEwt.,  for  all  J  and  proceeds  as  in  the  unswept  trailing  edge. 

MIN 


A. 3  ITERATION  PROCEDURE 

For  the  first  set  of  potential  flow  solution  iterations,  approximate 
wakes  composed  of  straight  line  vortices  trailing  downstream  from  the  canard 
and  wing  trailing  edges  are  constructed  using  the  slope  at  the  trailing  edge 
and  a  streamwise  decay  of  the  slope.  Subsequently,  after  each  set  of  itera¬ 
tions  the  wake  shape  is  updated  to  only  one  more  downstream  station  beyond  the 
previous  update,  i.e.,  after  the  latest  set  of  iterations  the  wake  is  "updat- 


wn  ~w  TTwrm  w  w 


ed"  from  the  trailing  edge  to  I  =*  IM.  Beyond  this  point  each  vortex  line  is 
continued  back  with  a  slope  at  each  downstream  segment  of  half  the  slope  of 
the  previous  segment.  This  scheme  of  slope  decay  prevents  the  wake  geometry 
from  diverging  or  becoming  chaotic  but  does  not  significantly  affect  the  solu¬ 
tion  in  the  upstream  regions  while  they  are  iteratively  building  up.  The 
above  procedure  is  illustrated  in  Figure  26.  Typically,  a  set  of  eight  itera¬ 
tions  is  performed  between  each  wake  update. 

A. A  WAKE  ROLL-UP 

It  is  well  known  that  trailing  vortex  sheets  eventually  roll  up  as  they 
proceed  downstream.  The  procedure  employed  here  for  tracking  vortex  lines 
allows  this  phenomenon.  As  a  result,  two  serious  difficulties  are  encountered 
in  the  computational  scheme. 

First,  the  flow  solution  uses  a  computational  space  in  which  the  non~ 

planar  vortex  sheets  of  the  canard  and  wing  are  mapped  onto  parallel  planar 

surfaces.  When  wake  roll-up  occurs,  the  spanwise  coordinate,  YV  ,  of  the 
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sheet  becomes  multiple  valued  and  the  mapping  function  breaks  down.  Second, 
the  discretization  of  the  wake  may  result  in  the  interval  between  vortex  lines 
becoming  large  compared  to  the  spanwise  radius  of  curvature  of  the  wake.  The 
"tracking"  procedure  then  becomes  invalid,  producing  a  chaotic  wake  surface. 

When  wake  roll-up  occurs,  the  vortex  sheet  is  replaced  by  one  which  has 
single  valued  coordinates  and  a  vorticity  distribution  which  approximates  that 
of  the  rolled  up  wake.  This  is  accomplished  by  monitoring  the  spanwise  sepa¬ 
ration  between  vortex  lines  at  each  streamwise  station  I  as  the  wake  is  being 

updated.  Roll-up  begins  as  values  of  YV  approach  each  other.  Each  down- 
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stream  station  I  searches  for  roll-up  with  a  "proximity  test": 

YVI  l+1  "  YVI  L-  DYTEST  (A-2) 

where  DYTEST  is  a  set  parameter  (1/2  of  the  first  spanwise  grid  interval). 

When  roll-up  is  detected,  the  two  vortex  lines  are  "merged"  to  a  single 
line  of  strength  dr^+dr^j  located  at  the  centroid  of  the  magnitudes  of  their 
strength: 
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Fijpir*  26.  Waka  shapa  ftarativa  devalopmant. 
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dr 
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ZVI,L.|  drL*l  +  ZVI,L  drL 
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Merged  ldrL  +|drL+i 


(A-3b) 


Scanning  over  L  at  a  given  station  I  is  repeated  until  the  wake  coordi¬ 
nates  are  "single  valued,"  and  the  vortex  lines  are  renumbered  to  have  consec¬ 
utive  indices.  The  vortex  line  tracking  (Equation  A-l)  then  continues  down¬ 
stream  as  before.  The  vortex  line  at  the  symmetry  plane  or  root  juncture  is 
continued  straight  back  and  never  "merged."  The  proximity  test  commences  with 
L  =  2.  The  merging  process  is  illustrated  for  a  single  surface  in  Figure  27. 


A. 5  TREATMENT  OF  LARGE  SPANWISE  WAKE  SLOPE 


Even  though  roll-up  is  prevented  from  occurring  it  is  possible  for  the 
spanwise  slope  of  the  wake  to  become  large  and  to  undergo  large  changes  in  the 
spanwise  direction.  Both  of  these  circumstances  have  been  found  to  induce 
numerical  instability  in  the  flow  solution.  It  then  becomes  desirable  to  sup¬ 
press  local  instances  of  such  occurrences  especially  in  the  far  wake  in  order 
to  preserve  the  solution  in  the  near  wake  region  where  it  is  of  importance. 
This  is  accomplished  by  calculating  the  spanwise  slope  of  the  wake  at  a  sta¬ 
tion  I  after  any  necessary  vortex  merging  has  been  carried  out  and  testing  the 
maximum  slope  against  a  reference  criterion  SLTEST  (set  at  0.50).  If  the  max¬ 
imum  slope  is  found  between  indices  J  =  JKM  -  1  and  J  “  JKM 


AZ 

ZI,JKM  ZI,JKM-1 

Ay 

Y  -  Y 

JKM  JKM-1 

>  SLTEST 


(A-Aa) 


a  reference  line  is  established 

ZREF  =  **  (ZI,JKM  +  ZI,JKM-1) 


(A~4b) 


and  the  wake  cross-section  is  scaled  about  this  line: 


NEW 


Ay 

SLTEST  E 


(A-4c) 


A  modified  Z___(Y)  is  introduced  in  the  case  of  swept  trailing  edges  at 

Rtr 

streamwise  stations  that  intercept  the  planform.  These  procedures  are  illus¬ 
trated  in  detail  in  Figure  28a  and  Figure  28b.  In  the  latter  case  ZREp(Y)  is 
chosen  to  insure  that  the  vortex  sheet  makes  contact  with  the  trailing  edge 
after  scaling.  If  contact  is  made  with  only  one  edge  as  for  the  inner  section 
in  Figure  28b,  Z  is  continued  horizontally  from  the  trailing  edge. 

Kbf 


A. 6  PREVENTION  OF  CROSSING  WAKES 

Fluid  dynamic  principles  rule  out  the  crossing  of  wing  and  canard  wakes 
except  along  lines  of  zero  vorticity  or  along  stagnation  lines.  However,  due 
to  finite  computational  grid  intervals  and  the  interpolation  procedures  em¬ 
ployed  in  the  numerical  solution,  the  wakes  may  be  found  to  cross  when  in  fact 
they  should  be  simply  in  close  proximity.  At  this  point  the  computational 
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procedure  can  be  expected  to  break  down.  To  prevent  this,  the  wake  is  tested 
by  scanning  across  spanwise  grid  points  J  at  a  fixed  streamwise  station  I  af¬ 
ter  vortex  merging  and  slope  corrections  have  been  made.  If  "crossing"  is 
detected  at  any  J  the  canard  wake  is  thereafter  continued  parallel  to  the  wing 
wake  with  the  separation  equal  to  that  at  station  1-1.  In  fact,  numerical 
difficulties  are  encountered  when  the  wakes  approach  each  other  too  closely 
before  crossing  actually  occurs.  Therefore,  the  criterion  used  to  detect 
"crossing"  is  to  compare  the  wake  separation  with  a  minimum  acceptable  sepa¬ 
ration  SFIX,  which  is  taken  to  be  0.20  of  the  original  wing-canard  vertical 
separation.  Once  the  canard  wake  is  continued  back  in  parallel  fashion,  the 
vortex  line  generators  are  assumed  to  move  straight  back  with  constant  span- 
wise  distribution  of  circulation  and  potential  jump. 

A. 7  WAKE  BOUNDARY  CONDITIONS 

Once  vortex  merging  and  slope  corrections  have  been  made,  the  new  distri¬ 
bution  of  circulation,  Tl^,  is  found  at  each  streamwise  station  for  all  points 
by  simple  numerical  integration  of  the  vortex  strength  dr^.  These  values  of 
Tl  are  interpolated  to  determine  at  the  vortex  line  location  and  further 
interpolated  to  determine  T  at  the  grid  locations.  The  wake  coordinates  at 

i  «  J 

the  grid  points  are  also  determined  by  interpolation  of  the  vortex  line 
coordinates.  The  boundary  condition  on  the  wing  and  canard  wake  surfaces  is 
the  jump  in  potential  equal  to  T  .  Beyond  the  last  wake  update  point  no 
change  is  considered  to  occur  in  the  spanwise  distribution  of  T.  To  prevent 
numerical  "kinks"  in  wake  shape  at  the  downstream  boundary,  the  slope  in  the  x 
direction  of  the  wake  is  held  constant  for  the  three  streamwise  stations  near¬ 
est  the  boundary. 

The  vertical  grid  construction  considers  the  wake  surfaces  to  extend  be¬ 
yond  the  wing  and  canard  tip  stations  to  the  outer  Y  boundary.  Thus,  beyond 
the  most  outboard  vortex  line  of  the  wake  the  grid  lines  are  continued  with 
constant  values  of  Z.  The  edge  of  the  wake  is  free  to  expand  or  shrink  as  the 
outermost  vortex  line  is  tracked  downstream. 


APPENDIX  B 

GLOSSARY  OF  OUTPUT  VARIABLES 


CCD/ CAV 

CCF/CAV 

CCL/CAV 

CCM/CAV/MAC 

CCMAX 

CCAV 

CD 

CDINT 

CF 

CIR 

CL 

CLINT 

CM 

CMLOC 

CP 

DRAG 

ETA 

I 

IL 

INOSE 

INOSEC 

IT 

ITAIL 


Spanload  drag. 

Spanload  skin  friction. 

Spanload  lift. 

Spanload  pitching  moment. 

Maximum  correction  to  flow  potential  for  the  current  itera¬ 
tion. 

Average  of  the  absolute  value  of  all  the  flowfield  correc¬ 
tions  for  the  current  iteration. 

Drag  coefficient  (local  or  total,  as  appropriate). 

Local  drag  coefficient  in  spanload  tables. 

Skin  friction  coefficient. 

Circulation. 

Lift  coefficient  (local  or  total  as  appropriate). 

Local  lift  coefficient  in  spanload  tables. 

Pitching  moment  coefficient. 

Local  pitching  moment  coefficient  (about  quarter  chord 
point) . 

Pressure  coefficient. 

Local  drag  contribution  in  the  body  force  output. 

Spanwise  computational  ordinate. 

General  streamwise  index. 

Streamwise  index  of  leading  edge  point  in  crude  grid. 

Streamwise  index  of  the  first  point  on  the  body  in  the  body 
fine  grid. 

Streamwise  index  of  the  first  point  on  the  body  in  the  crude 
grid. 

Streamwise  index  of  the  trailing  edge  point  at  an  analysis 
station  in  the  crude  grid. 

Streamwise  index  of  the  last  point  on  the  body  in  the  body 
fine  grid. 


ITAILC 

ITER 

J 

JSD 

KLO 

LOAD 

NSP 

PCT 

PHI 

RMAX 

RSD 

RSDAV 

SEXPC 

SHIFT 

U 

V 
W 

WCORD 

X 

XI 
XLE 
XNOSE 
XTAIL 
XTE 
XWF 
X/C 

Y 

Y/C 

Z 

ZETA 

2Y/B 


Streamwise  index  of  the  last  point  on  the  body  in  the  crude 
grid. 

Iteration  count. 

General  spanwise  index. 

Spanwise  index  of  first  solution  plane  outside  the  body  com¬ 
putational  surface. 

Vertical  index  for  bottom  of  body  computational  surface. 

Local  lift  contribution  in  body  force  output. 

Number  of  supersonic  points. 

Fractional  cord  distance  from  leading  or  trailing  edge  during 
crude  grid  shift  search. 

Perturbation  velocity  potential. 

Maximum  body  radius. 

Residual  with  maximum  absolute  value. 

Average  of  the  absolute  value  of  the  residual  at  all  grid 
points  for  the  current  iteration. 

Canard  exposed  area. 

Grid  shift  parameter. 

Streamwise  perturbation  velocity. 

Spanwise  perturbation  velocity. 

Vertical  perturbation  velocity  or  over-relaxation  factor. 
Local  cord  for  wing  or  canard. 

Streamwise  physical  coordinate. 

Streamwise  computational  coordinate. 

Wing  or  canard  leading  edge  X-location. 

Body  nose  X-location. 

Body  tail  X-location. 

Wing  or  canard  trailing  edge  location. 

Wing  or  canard  fine  grid  streamwise  physical  coordinate. 
Percent  fraction  of  local  chord. 

Spanwise  physical  coordinate. 

Percent  chord  wing  or  canard  section  ordinate. 

Vertical  physical  coordinate. 

Vertical  computational  coordinate. 

Fraction  of  semispan  (on  wing  or  canard  as  appropriate). 
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